{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from [this Jupyter notebook](https://nbviewer.jupyter.org/github/krischer/seismo_live/blob/master/notebooks/Computational%20Seismology/The Finite-Difference Method/fd_taylor_operators_advanced.ipynb) by Heiner Igel ([@heinerigel](https://github.com/heinerigel)), Lion Krischer ([@krischer](https://github.com/krischer)) which is a supplemenatry material to the book [Computational Seismology: A Practical Introduction](http://www.computational-seismology.org/),  additional modifications by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Improving the accuracy of FD operators based on Taylor series expansion\n",
    "\n",
    "After estimating the accuracy of FD operators by calculating the remainders of Taylor series expansions, the next obvious question is how could we improve the accuracy of FD approximations.\n",
    "\n",
    "## Estimation of FD operators by Taylor series expansion\n",
    "\n",
    "In [this lession](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/03_Intro_finite_differences/1_fd_intro.ipynb), we introduced the 3-point FD approximation to the second derivative of a function:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx \\frac{f(x+dx)-2f(x)+f(x-dx)}{dx^2} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "from a Taylor series expansion. Here, I want to introduce a more systematic approach to derive this operator, which is required for a later improvement of the FD operator accuracy.\n",
    "\n",
    "The Taylor expansion of $f(x + dx)$ around $x$ is defined as \n",
    "\n",
    "$$\n",
    "f(x+dx)=\\sum_{n=0}^\\infty \\frac{f^{(n)}(x)}{n!}dx^{n}  \n",
    "$$\n",
    "\n",
    "Finite-difference operators can be calculated by seeking weights (here: $a$, $b$, $c$) with which function values have to be multiplied to obtain an interpolation or a derivative. \n",
    "\n",
    "### Example: second derivative\n",
    "\n",
    "We want to approximate the second derivative of a function based on function values at the points $f(x+dx)$, $f(x)$ and $f(x-dx)$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx a ~ f(x+dx) + b ~ f(x) + c ~ f(x-dx)\n",
    "\\end{equation}\n",
    "\n",
    "The individual terms in the derivative can be expressed by the following Taylor series expansion up to the quadratic term:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "a ~ f(x + dx) & \\ = \\ a ~ \\left[ ~  f(x) +  \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{1}{2!}  \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2   + \\dotsc  ~ \\right] \\nonumber\\\\\n",
    "b ~ f(x) & \\  = \\ b ~ \\left[ ~  f(x)  ~ \\right] \\nonumber\\\\\n",
    "c ~ f(x - dx) & \\ = \\ c ~ \\left[ ~  f(x) -  \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{1}{2!}  \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2   - \\dotsc  ~ \\right] \\nonumber\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Inserting these approximations into eq. (1) and sorting of the terms leads to \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx a ~ f(x+dx) + b ~ f(x) + c ~ f(x-dx) \\approx \\nonumber\n",
    "\\end{equation}\n",
    "$$\n",
    "\\begin{align}\n",
    "f(x) ~ \\biggl[&a  ~~+~~            ~~~~b           &+~~  c\\biggr] \\notag\\\\\n",
    "+dx ~ \\frac{\\partial f(x)}{\\partial x} ~ \\biggl[&a  ~~\\phantom{+}~~ \\phantom{b}  &-~~  c\\biggr] \\notag\\\\\n",
    "+\\frac{dx^2}{2!} ~ \\frac{\\partial^2 f(x)}{\\partial x^2} ~ \\biggl[&a  ~~\\phantom{+}~~ \\phantom{b}  &+~~  c\\biggr] \\notag\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The RHS of this equation can only be an approximation to $\\frac{\\partial^2 f(x)}{\\partial x^2}$, if the following conditions are fulfilled:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&a  ~~+~~            ~~~~b           &+~~  c & = & 0 \\notag\\\\\n",
    "&a  ~~\\phantom{+}~~ \\phantom{b}  &-~~  c & = & 0 \\notag\\\\\n",
    "&a  ~~\\phantom{+}~~ \\phantom{b}  &+~~  c & = & \\frac{2!}{\\mathrm{d}x^2}\\notag\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This system of linear equations can be expressed in matrix form as:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "  1 & 1 & 1 \\\\\n",
    "  1 & 0 & -1 \\\\\n",
    "  1 & 0 & 1\n",
    " \\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "  a\\\\\n",
    " b \\\\\n",
    " c\n",
    " \\end{pmatrix}\n",
    " =\\begin{pmatrix}\n",
    "  0\\\\\n",
    " 0 \\\\\n",
    " \\frac{2!}{dx^2}\n",
    " \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "or \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{Ax ~= ~s}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "with the system matrix $\\mathbf{A}$, the vector $\\mathbf{x}$ containing the desired FD weights $a,\\; b,\\; c$ and the solution vector $\\mathbf{s}$. Formally, the FD operator weights can be estimated by matrix inversion:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{x ~= A^{-1}~s}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Let's try to estimate the FD weighting coefficients by symbolic matrix inversion with `SymPy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import SymPy libraries\n",
    "from sympy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(dx**(-2), -2/dx**2, dx**(-2))}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define symbols\n",
    "a, b, c, dx = symbols('a b c dx')\n",
    "\n",
    "# In SymPy we define the matrix equation by a matrix consisting of the \n",
    "# system matrix A and an additional column defining the solution vector s\n",
    "system = Matrix(([1, 1, 1, 0], [1, 0, -1, 0], [1, 0, 1, 2/dx**2]))\n",
    "\n",
    "# Symbolic solution of the linear system can be easily obtained by linsolve\n",
    "linsolve(system, (a, b, c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the resulting weighting coefficients are:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "a \\\\\n",
    "b\\\\\n",
    "c\n",
    "\\end{pmatrix}\n",
    "=\\begin{pmatrix}\n",
    "\\frac{1}{\\mathrm{d}x^2} \\\\\n",
    "-\\frac{2}{\\mathrm{d}x^2} \\\\\n",
    "\\frac{1}{\\mathrm{d}x^2}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Inserting the coefficients $a,\\; b,\\; c$ in eq. (1) leads to \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx \\frac{f(x+dx)-2f(x)+f(x-dx)}{dx^2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "which is the well known 3-point operator for the 2nd derivative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improving the accuracy of the 2nd derivative FD approximation\n",
    "\n",
    "We can use the above approach to estimate higher order and therefore more accurate FD operators. Let's assume we want to approximate the second derivative by a 5- instead of a 3-point operator, at the points $f(x+2dx)$, $f(x+dx)$, $f(x)$, $f(x-dx)$ and $f(x-2dx)$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx a ~ f(x+2dx) + b ~ f(x+dx) + c ~ f(x) + d ~ f(x-dx) + e ~ f(x-2dx) \\notag\n",
    "\\end{equation}\n",
    "\n",
    "To solve for the 5 unknown FD weights, we have to expand 5 terms of the Taylor series, which leads to a higher error order and therefore a more accurate FD operator compared to the 3-point operator. The resulting Taylor series expansions are:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "a ~ f(x + 2dx) & \\ = \\ a ~ \\left[ ~  f(x) +  2 \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{4}{2}  \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2 +  \\frac{8}{6}  \\frac{\\partial^3 f(x)}{\\partial x^3} dx^3 +  \\frac{16}{24}  \\frac{\\partial^4 f(x)}{\\partial x^4} dx^4   + \\dotsc  ~ \\right] \\nonumber\\\\\n",
    "b ~ f(x + dx) & \\ = \\ b ~ \\left[ ~  f(x) +  \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{1}{2} \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2 +  \\frac{1}{6} \\frac{\\partial^3 f(x)}{\\partial x^3} dx^3 +  \\frac{1}{24}  \\frac{\\partial^4 f(x)}{\\partial x^4} dx^4   + \\dotsc  ~ \\right] \\nonumber\\\\\n",
    "c ~ f(x) & \\  = \\ c ~ \\left[ ~  f(x)  ~ \\right] \\nonumber\\\\\n",
    "d ~ f(x - dx) & \\ = \\ d ~ \\left[ ~  f(x) -  \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{1}{2} \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2 -  \\frac{1}{6} \\frac{\\partial^3 f(x)}{\\partial x^3} dx^3 +  \\frac{1}{24}  \\frac{\\partial^4 f(x)}{\\partial x^4} dx^4   + \\dotsc  ~ \\right] \\nonumber\\\\\n",
    "e ~ f(x - 2dx) & \\ = \\ e ~ \\left[ ~  f(x) -  2 \\frac{\\partial f(x)}{\\partial x} dx +  \\frac{4}{2}  \\frac{\\partial^2 f(x)}{\\partial x^2} dx^2 -  \\frac{8}{6} \\frac{\\partial^3 f(x)}{\\partial x^3} dx^3 +  \\frac{16}{24} \\frac{\\partial^4 f(x)}{\\partial x^4} dx^4   + \\dotsc  ~ \\right] \\nonumber\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Similar to the 3-point operator, we insert the terms into the 5-point FD approximation and sort the terms\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(x)}{\\partial x^2} \\approx a ~ f(x+2dx) + b ~ f(x+dx) + c ~ f(x) + d ~ f(x-dx) + e ~ f(x-2dx)  \\approx \\nonumber\n",
    "\\end{equation}\n",
    "$$\n",
    "\\begin{align}\n",
    "f(x) ~ \\biggl[&a  ~~+~~   b   ~~+~~ c ~~+~~  d ~~+~~  e~~           \\biggr] \\notag\\\\\n",
    "+dx ~ \\frac{\\partial f(x)}{\\partial x} ~ \\biggl[&2a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~-~~  d ~~-~~  2e~~           \\biggr] \\notag\\\\\n",
    "+\\frac{dx^2}{2} ~ \\frac{\\partial^2 f(x)}{\\partial x^2} ~ \\biggl[&4a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~+~~  d ~~+~~  4e~~           \\biggr] \\notag\\\\\n",
    "+\\frac{dx^3}{6} ~ \\frac{\\partial^3 f(x)}{\\partial x^3} ~ \\biggl[&8a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~-~~  d ~~-~~  8e~~           \\biggr] \\notag\\\\\n",
    "+\\frac{dx^4}{24} ~ \\frac{\\partial^4 f(x)}{\\partial x^4} ~ \\biggl[&16a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~+~~  d ~~+~~  16e~~           \\biggr] \\notag\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "leading to the following conditions:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\biggl[&a  ~~+~~   b   ~~+~~ c ~~+~~  d ~~+~~  e~~           \\biggr] ~=~ 0 \\notag\\\\\n",
    "\\biggl[&2a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~-~~  d ~~-~~  2e~~           \\biggr] ~=~ 0 \\notag\\\\\n",
    "\\biggl[&4a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~+~~  d ~~+~~  4e~~           \\biggr] ~=~ \\frac{2}{dx^2} \\notag\\\\\n",
    "\\biggl[&8a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~-~~  d ~~-~~  8e~~           \\biggr] ~=~ 0\\notag\\\\\n",
    "\\biggl[&16a  ~~+~~   b   ~~\\phantom{+}~~ \\phantom{c} ~~+~~  d ~~+~~  16e~~           \\biggr] ~=~ 0\\notag\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "or as matrix equation:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "  1 & 1 & 1 & 1 & 1 \\\\\n",
    "  2 & 1 & 0 & -1 & -2 \\\\\n",
    "  4 & 1 & 0 & 1 & 4 \\\\\n",
    "  8 & 1 & 0 & -1 & -8 \\\\\n",
    "  16 & 1 & 0 & 1 & 16 \\\\\n",
    " \\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "  a\\\\\n",
    " b \\\\\n",
    " c \\\\\n",
    " d \\\\\n",
    " e \\\\\n",
    " \\end{pmatrix}\n",
    " =\\begin{pmatrix}\n",
    "  0\\\\\n",
    " 0 \\\\\n",
    " \\frac{2!}{dx^2} \\\\\n",
    " 0 \\\\\n",
    " 0\n",
    " \\end{pmatrix}\n",
    " \\nonumber\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(-1/(12*dx**2), 4/(3*dx**2), -5/(2*dx**2), 4/(3*dx**2), -1/(12*dx**2))}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define symbols\n",
    "a, b, c, d, e, dx = symbols('a b c d e dx')\n",
    "\n",
    "# In SymPy we define the matrix equation by a matrix consisting of the \n",
    "# system matrix A and an additional column defining the solution vector s\n",
    "system = Matrix(([1, 1, 1, 1, 1, 0], [2, 1, 0, -1, -2, 0], [4, 1, 0, 1, 4, 2/dx**2], [8, 1, 0, -1, -8, 0], [16, 1, 0, 1, 16, 0]))\n",
    "\n",
    "# Symbolic solution of the linear system can be easily obtained by linsolve\n",
    "linsolve(system, (a, b, c, d, e))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the resulting weighting coefficients for the 5-point FD operator are:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "a \\\\\n",
    "b\\\\\n",
    "c\\\\\n",
    "d\\\\\n",
    "e\n",
    "\\end{pmatrix}\n",
    "=\\begin{pmatrix}\n",
    "-\\frac{1}{12\\mathrm{d}x^2} \\\\\n",
    " \\frac{4}{3\\mathrm{d}x^2} \\\\\n",
    "-\\frac{5}{2\\mathrm{d}x^2}\\\\\n",
    "\\frac{4}{3\\mathrm{d}x^2}\\\\\n",
    "-\\frac{1}{12\\mathrm{d}x^2}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "##### Exercise\n",
    "\n",
    "To investigate the effect of the 5-point operator on the finite-difference based solution of the 1D acoustic wave equation:\n",
    "\n",
    "* Implement the 5-point operator to approximate the spatial 2nd derivative in the 1D acoustic FD modelling code. Add an option to use the 3-point operator\n",
    "* Compare the impact of the 3-point and 5-point operators on the modelling result by plotting the analytical seismograms $d^{analyt}$ and FD seismograms $d^{FD}$ as well as their differences, respectively. Use $v_0 = 333\\; m/s$, $dx = 2\\; m$ and $dt = 1\\;ms$ as modelling parameters. Quantify the data misfit also in terms of the L2-norm of the data residuals:\n",
    "\\begin{equation}\n",
    "E = \\frac{1}{2}\\sum_{i=0}^{nt}\\biggl(d^{FD}_i-d^{analyt}_i\\biggr)^2\\nonumber\n",
    "\\end{equation}\n",
    "where $nt$ denotes the number of time samples\n",
    "* Benchmark and compare the runtimes of the 3-point and 5-point operator by using the magic function `%timeit` and the same spatio-temporal model discretizations $dx$ and $dt$ in order to achieve comparable results\n",
    "* Estimate the factor $\\zeta$ in the Courant criterion:\n",
    "\\begin{equation}\n",
    "dt \\le \\frac{dx}{\\zeta v_{max}}\\nonumber\n",
    "\\end{equation}\n",
    "for the 5-point operator by forward modelling. Start with the modelling parameters $v_{max} = 333\\; m/s$ and $dx = 0.5\\; m$ and the Courant criterion for the 3pt operator ($\\zeta=1$). Vary $\\zeta$ during different forward modelling runs, until you resolved the transition between stable and unstable FD solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "xmax = 500 # maximum spatial extension of the 1D model (m)\n",
    "dx   = 0.5 # grid point distance in x-direction\n",
    "\n",
    "tmax = 1.001   # maximum recording time of the seismogram (s)\n",
    "dt   = 0.0010  # time step\n",
    "\n",
    "vp0  = 333.   # P-wave speed in medium (m/s)\n",
    "\n",
    "# acquisition geometry\n",
    "xr = 365.0 # receiver position (m)\n",
    "xsrc = 249.5 # source position (m)\n",
    "\n",
    "f0   = 25. # dominant frequency of the source (Hz)\n",
    "t0   = 4. / f0 # source time shift (s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1D Wave Propagation (Finite Difference Solution) \n",
    "# ------------------------------------------------\n",
    "def FD_1D_acoustic(dt,dx,op):\n",
    "        \n",
    "    nx = (int)(xmax/dx) # number of grid points in x-direction        \n",
    "    nt = (int)(tmax/dt) # maximum number of time steps            \n",
    "    \n",
    "    ir = (int)(xr/dx)      # receiver location in grid in x-direction    \n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "\n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of a Gaussian\n",
    "    src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "    # Analytical solution\n",
    "    # -------------------\n",
    "    G    = time * 0.\n",
    "\n",
    "    # Initialize coordinates\n",
    "    # ----------------------\n",
    "    x    = np.arange(nx)\n",
    "    x    = x * dx       # coordinate in x-direction\n",
    "\n",
    "    for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "        if (time[it] - np.abs(x[ir] - x[isrc]) / vp0) >= 0:\n",
    "            G[it] = 1. / (2 * vp0)\n",
    "    seis_ana   = np.convolve(G, src * dt)\n",
    "    seis_ana   = seis_ana[0:nt]\n",
    "    \n",
    "    # Initialize empty pressure arrays\n",
    "    # --------------------------------\n",
    "    p    = np.zeros(nx) # p at time n (now)\n",
    "    pold = np.zeros(nx) # p at time n-1 (past)\n",
    "    pnew = np.zeros(nx) # p at time n+1 (present)\n",
    "    d2px = np.zeros(nx) # 2nd space derivative of p\n",
    "\n",
    "    # Initialize model (assume homogeneous model)\n",
    "    # -------------------------------------------\n",
    "    vp    = np.zeros(nx)\n",
    "    vp    = vp + vp0       # initialize wave velocity in model\n",
    "\n",
    "    # Initialize empty seismogram\n",
    "    # ---------------------------\n",
    "    seis = np.zeros(nt)    \n",
    "    \n",
    "    # Calculate Partial Derivatives\n",
    "    # -----------------------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # FD approximation of spatial derivative by 3 point operator\n",
    "        if(op==3):\n",
    "            for i in range(1, nx - 1):\n",
    "                d2px[i] = (p[i + 1] - 2 * p[i] + p[i - 1]) / dx ** 2\n",
    "        \n",
    "        # ADD YOUR 5-POINT FD OPERATOR HERE \n",
    "        #if(op==5):\n",
    "                \n",
    "        # Time Extrapolation\n",
    "        # ------------------\n",
    "        pnew = 2 * p - pold + vp ** 2 * dt ** 2 * d2px\n",
    "\n",
    "        # Add Source Term at isrc\n",
    "        # -----------------------\n",
    "        # Absolute pressure w.r.t analytical solution\n",
    "        pnew[isrc] = pnew[isrc] + src[it] / dx * dt ** 2\n",
    "                \n",
    "        # Remap Time Levels\n",
    "        # -----------------\n",
    "        pold, p = p, pnew\n",
    "    \n",
    "        # Output of Seismogram\n",
    "        # --------------------\n",
    "        seis[it] = p[ir]\n",
    "        \n",
    "    return time, seis, seis_ana"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = 2.0 # grid point distance in x-direction (m)\n",
    "dt = 0.0010  # time step (s)\n",
    "\n",
    "# run FD code with 3-point operator\n",
    "op = 3 # use 3-point operator\n",
    "time, seis_3pt, seis_ana = FD_1D_acoustic(dt,dx,op)\n",
    "\n",
    "# run FD code with 5-point operator\n",
    "op = 5 # use 5-point operator\n",
    "# RUN FD CODE WITH 5-POINT OPERATOR HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvoAAAFNCAYAAABxOMuFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl4VdXZ9/HvTRgSBgkFFAQEFAhkIGEWZVKqoiKoBUGtAmoRlbY+VRza51WKQ9U6UB+tU1W0VUBRKSJ1QIiIghIgQCDMBqGgIAgEkkCG9f5xdg6H5GSADCckv891nSt7Wmvd+2QT7rPPWmubcw4REREREaleaoU6ABERERERKX9K9EVEREREqiEl+iIiIiIi1ZASfRERERGRakiJvoiIiIhINaREX0RERESkGlKiLyISYmZ2yMzODnUccmLMbK2ZDQp1HCIiRTHNoy8iUnZm1g94AogBcoFU4E7n3LKQBiYiIjVW7VAHICJyqjOz04C5wG3AO0BdoD9wJJRxlRczC3PO5YY6jqrCzGo753JCHYeISEnUdUdEpOw6ATjnpjvncp1zmc65T51zq/MPMLObzCzVzH42s0/MrG3APmdmHbzly8xsnZmlm9l/zexub/sgM9thZveY2W4z22VmV3rHbzSzfWb2x4A665nZVDPb6b2mmlm9gP33eHXsNLNbCsQwzcxeMLN5ZnYYuMDMLjezlWZ20My2m9nkgLraeeXHeft+NrMJZtbLzFab2X4ze66oN8/Mwszsj2a2xTvv5WbWxtt3npktM7MD3s/zAsolmtnDZva11/3pQzNramZveXEuM7N2Bd7n35nZVjP7ycz+ama1vH3nmNkCM9vr7XvLzCIDyqaZ2b1mtho4bGa1vW2/9Pb3NrMkr90fzezpgLLDvG4++72YuxSo927vfTpgZjPNLLzYq01EpJSU6IuIlN1GINfM3jCzS82sSeBOM7sS+CNwNdAc+BKYXkRdrwK3OucaAbHAgoB9LYBwoBXwAPAK8GugB75vEB4I6Ov/J+BcIAGIB3oD/+vFMwT4A/BLoAMwMEgc1wGPAI2AxcBh4EYgErgcuM07r0B9gI7AKGCqF8Mv8XVnusbMgrWDF8u1wGXAacBNQIaZ/QL4CHgWaAo8DXxkZk0Dyo4GbvDek3OAJcDrwC/wdZ96sEBbVwE9ge7AcK8tAAP+ApwJdAHaAJMLlL3WO/fIIHf0/wb8zTl3mhfHOwBm1gnf7/pOfL/7ecCHZlY3oOw1wBCgPdAVGFvE+yQickKU6IuIlJFz7iDQD3D4ku89ZjbHzM7wDrkV+ItzLtVLEB8FEgLv6gfIBqLN7DTn3M/OuRUF9j3inMsGZgDN8CWX6c65tcBafIkiwPXAFOfcbufcHuDP+BJi8CWWrzvn1jrnMrx9Bf3bOfeVcy7POZflnEt0zq3x1lfjS14LJu4Pecd+iu+DwXSv/f/i+3DTrYi38Bbgf51zG5zPKufcXnxJ9Sbn3D+dcznOuenAeuCKgLKvO+e2OOcOAP8Btjjn5nvv87tB2nzcObfPOfc9vg8j1wI45zY75z5zzh3x3q+ng5zfs8657c65zCDnkA10MLNmzrlDzrml3vZRwEde3dnAk0AEcF5A2Wedczudc/uAD/F9OBMRKTMl+iIi5cBL4sc651rjuxN/Jr5EEqAt8Dev68Z+YB++O8itglT1K3x3treZ2Rdm1jdg396AvvL5yeaPAfszgYbe8pnAtoB927xt+fu2B+wLXA66zcz6mNlCM9tjZgeACfg+aAQqGEtRsRXUBtgSZHvBc8BbD3zfTrTNwPPyvydmdrqZzTBfd6mDwL8ofH7B3qd8N+PrwrXe6zI0NNg5OOfyvHoCz+GHgOWMIDGLiJwUJfoiIuXMObcemIYv4QdfYnercy4y4BXhnPs6SNllzrnhwOnAbLwuICdhJ74PGPnO8rYB7AJaB+xrE+w0Cqy/DcwB2jjnGgMv4vuwUh624+vuUlDBcwDfefy3DG0Fnmvge/IXfOfc1et+82sKn1+R09Q55zY5567F93t7HJhlZg0ocA5mZl4MZTkHEZFSUaIvIlJGZtbZzO4ys9beeht8XULyu2+8CNxvZjHe/sZmNjJIPXXN7Hoza+x18ziIb6rOkzEd+F8za25mzfD16f+Xt+8dYJyZdTGz+t6+kjQC9jnnssysN74+/OXlH8BDZtbRfLp6/fDnAZ3M7Dpv8OsoIBrfDEcna5KZNfF+R78HZnrbGwGHgP1m1gqYdCKVmtmvzay5d8d+v7c5F997fbmZDTazOsBd+GZjKvQhT0SkvCnRFxEpu3R8A1G/Md8sNUuBFHxJHc65D/Dd5Z3hdQtJAS4toq4bgDTvuAn47iyfjIeBJGA1sAZY4W3DOfcffANcFwKb8Q1gheKnA70dmGJm6fg+GJzsNw3BPO3V9ym+DzevAhFeP/2h+N7HvcA9wFDn3E9laOvfwHIgGd9A31e97X/GN0D3gLf9/ROsdwiw1swO4RuYO9obr7AB3+/w/4Cf8I0vuMI5d7QM5yAiUip6YJaISA3nTfeYAtSrzvPDm5kDOjrnNoc6FhGRyqA7+iIiNZCZXeV1FWqC79uGD6tzki8iUhMp0RcRqZluBfbgm+0mF99TfUVEpBpR1x0RERERkWpId/RFRERERKohJfoiIiIiItVQ7VAHcKqKjIx0HTp0CHUYUoUcPnyYBg0ahDoMqWJ0XUgwui4kGF0XEszy5ct/cs41P5mySvRP0hlnnEFSUlKow5AqJDExkUGDBoU6DKlidF1IMLouJBhdFxKMmW072bLquiMiIiIiUg0p0RcRERERqYaU6IuIiIiIVEPqoy8iIiI1UnZ2Njt27CArKyvUoQDQuHFjUlNTQx2GhEh4eDitW7emTp065VanEn0RERGpkXbs2EGjRo1o164dZhbqcEhPT6dRo0ahDkNCwDnH3r172bFjB+3bty+3etV1R0RERGqkrKwsmjZtWiWSfKnZzIymTZuW+7dLSvRFRESkxlKSL1VFRVyLIU30zWyImW0ws81mdl+Q/fXMbKa3/xszaxew735v+wYzuyRg+2tmttvMUgrUNdnM/mtmyd7rspLqEhEREalIYWFhJCQkkJCQwPnnn09aWhqJiYk0btyYbt26ERUVxYABA5g7d26Z20pLSyM2NrbE4x599NHj1s8777wyt53vzjvvZNGiRQDcfPPNxMfH07VrV0aMGMGhQ4eKLZucnMy8efPKLZaymjZtGjt37jzhcnfffTcLFiyogIgKC1mib2ZhwPPApUA0cK2ZRRc47GbgZ+dcB+AZ4HGvbDQwGogBhgB/9+oDmOZtC+YZ51yC95pXirpERKQCOQe7dsGBA6GORCQ0IiIiSE5OJjk5ma+++op27doB0L9/f1auXMmGDRt49tlnmThxIp9//nmlxFQw0f/666/Lpd59+/axdOlSBgwYAMAzzzzDqlWrWL16NWeddRbPPfdcseVDkejn5uYWue9kEv3c3Fx++9vf8thjj5U1tFIJ5R393sBm59xW59xRYAYwvMAxw4E3vOVZwGDzfa8xHJjhnDvinPsO2OzVh3NuEbDvBOIosi4REak48+dDdDSceSY0aQLDhsFJ3BwTqfYSEhJ44IEHgibCX3zxhf8bgW7dupGeno5zjkmTJhEbG0tcXBwzZ84sVG7atGlMnDjRvz506FASExO57777yMzMJCEhgeuvvx6Ahg0bAhRZb/4TfUeMGEHnzp25/vrrcc4VanPWrFkMGXLsXuxpp53mrzczM9PfdWXs2LFMmDCB/v3706lTJ+bOncvRo0d54IEHmDlzJgkJCYXOKSsri3HjxhEXF0e3bt1YuHCh/zyHDx/OkCFDiIqK4s9//rO/zL/+9S969+5NQkICt956qz+pb9iwIQ888AB9+vRhyZIlTJkyhV69ehEbG8v48eNxzjFr1iySkpK4/vrrSUhIIDMzk88//5xu3boRFxfHTTfdxJEjRwBo164dU6ZMoV+/frz77ru0bduWvXv38sMPPxT7ey8PoUz0WwHbA9Z3eNuCHuOcywEOAE1LWTaYiWa22uve0+QE4hARkfKSlcXsWTkMGQLr1/s2OQcffgh9+yrZl5olP6lOSEjguuuuK/K47t27sz7/H0yAJ598kueff57k5GS+/PJLIiIieP/990lOTmbVqlXMnz+fSZMmsWvXrlLF89hjj/m/ZXjrrbeO21dcvStXrmTq1KmsW7eOrVu38tVXXxWq+6uvvqJHjx7HbRs3bhwtWrRg/fr1/Pa3v/VvT0tL44svvuCjjz5iwoQJ5OXlMWXKFEaNGkVycjKjRo06rp7nn38egDVr1jB9+nTGjBnjH9j67bff8tZbb5GcnMy7775LUlISqampzJw5k6+++ork5GTCwsL853v48GFiY2P55ptv6NevHxMnTmTZsmWkpKSQmZnJ3LlzGTFiBD179vTXa2aMHTuWmTNnsmbNGnJycnjhhRf88YWHh7N48WJGjx7t/30Ge4/KWyin1ww24qDgx7+ijilN2YJeAB7yjnsIeAq46UTqMrPxwHiA5s2bk5iYWEKTUpMcOnRI14QUouuigLw8zpk0mbzk06mV9za51KVWLUdengGOvd8f5tJLc3jmmWRqVePpInRdVA2NGzcmPT0dgNNOq7hpLQ8eTC9yX0REBF9++SXg69aRnp5ORkYGOTk5/tjAd83k5eUdtw2gZ8+e/P73v+eaa65h2LBhtGrVigULFnDVVVeRkZFB/fr1Oe+881i0aBExMTH+OrKysjh69Ki/vpycHDIyMvzrBdtJT08vst5GjRrRo0cPGjduzOHDh4mJiSE1NZX4+Pjj6ti+fTv169c/ru5nn32WZ555hrvvvps33niDX//612RnZzNs2DAOHz5MixYtaNu2LcuXLy8Uc6DExERuvfVW0tPTadWqFa1bt2blypVkZWUxaNAg6tatS05ODpdffjnz58+ndu3aJCUl+T94ZGZm+q+HsLAwLr74Yn878+bNY+rUqWRmZvLzzz/ToUMHBg0aRG5uLocPHyY9PZ01a9Zw1lln0bJlS9LT0xk5ciSvvPIKN998M845Lr/88uPijoyMZOvWrYXOJSsrq1z/NoQy0d8BtAlYbw0UvI+Tf8wOM6sNNMbXLac0ZY/jnPsxf9nMXgHyR7WUui7n3MvAywBRUVFu0KBBxTUpNUz+V5cigXRdFPDMM7DiS9oAsxjBPee8z4JFtUlOyuG/V95OC7eLYavnkJY2iJtuCnWwFUfXRdWQmppaKfPWl9RG/v78efTr169P7dq1jyu3ceNGYmJiCtX14IMPcvXVVzNv3jx++ctfMn/+fOrUqUN4eLj/2Dp16hAREUHDhg2pVasWjRo1omHDhse1kZOTQ/369f3rBdtp1KhRkfXWr1//uLLh4eHUqVMnaB1hYWFB348bbriBv/71r9x2223+evOPCwsLo2HDhoSHh1O3bt2g5cPCwo6LISwsjAYNGhQqU69ePSIiIqhVqxZjx47lL3/5S6G6wsPDiYyMBHyJ91133UVSUhJt2rRh8uTJOOf859KgQQP/7yzw3AJ/h2bGGWeccVzceXl5NGnSpNC5hIeH061bt0IxnaxQ3i9ZBnQ0s/ZmVhffgNg5BY6ZA4zxlkcAC5yv09ccYLQ3K097oCPwbXGNmVnLgNWrgPxZeU64LhEROQk//EDunx7wr26mA2++XZszW+Rx2UvD+Y17hSuYyyV8wiOPQE5OCGMVqUJWr17NQw89xB133FFo35YtW4iLi+Pee++lZ8+erF+/ngEDBjBz5kxyc3PZs2cPixYtonfv44cftmvXjuTkZPLy8ti+fTvffnss9alTpw7Z2dmF2ipNvcXp0qULmzdvBnz98gOXP/zwQzp37uw/9t133yUvL48tW7awdetWoqKiaNSoUdC7+fmx5Xe92bhxI99//z1RUVEAfPbZZ+zbt4/MzExmz57N+eefz+DBg5k1axa7d+8GfAOFt23bVqje/O4/zZo149ChQ8yaNcu/LzCezp07k5aW5j+nf/7znwwcOLDI92Ljxo2lmgGprEJ2R985l2NmE4FPgDDgNefcWjObAiQ55+YArwL/NLPN+O7kj/bKrjWzd4B1QA5wh3MuF8DMpgODgGZmtgN40Dn3KvCEmSXg65aTBtxaUl0iIlKOnnqKsEzf9HkpxLDimsf5Q2+AWtC6tf+wB/kz5229hPfeMwp0wxWpMEHGjobUl19+Sbdu3cjIyOD000/n2WefZfDgwYWOmzp1KgsXLiQsLIzo6GguvfRS6taty5IlS4iPj8fMeOKJJ2jRogVpaWn+cueffz7t27cnLi6O2NhYunfv7t83fvx4unbtSvfu3Y/rp3/VVVcFrTfY2IFgLr/8cl566SVuueUWnHOMGTOGgwcP4pwjPj7+uD7tUVFRDBw4kB9//JEXX3yR8PBwLrjgAh577DESEhK4//77j+unf/vttzNhwgTi4uKoXbs206ZNo169egD069ePG264gc2bN3PdddfRs2dPAB5++GEuvvhi8vLyqFOnDs8//zxt27Y9LubIyEh+85vfEBcXR7t27ejVq5d/X/6g4YiICJYsWcLrr7/OyJEjycnJoVevXkyYMCHo+5Cdnc3mzZv9cVQkCzYqWkoWFRXlNmzYEOowpArRV/ESjK4LT3o6uWe2JuzQQQCutH/z5MZhdOjg7d+xA845B44eBeAiPuXnHheRlBSieCuYrouqITU1lS5duoQ6DL/8rjvVWb9+/Zg7d66/a0wwY8eOZejQoYwYMaLM7U2bNo2kpKQSp+6sTB988AErVqzgoYceKrQv2DVpZsudcyf1qaAaD3USEZEq4/XX/Un+BjqRe+nQY0k++O7o33KLf3Uiz7F8OaxdW8lxikiFeuqpp/j+++9DHUZI5eTkcNddd1VKW0r0RUSkYjmHe/El/+pU7mT8hCD//fzP//gXL2MezdjD229XRoAiUln69OlD165diz1m2rRp5XI3H3zfDlSlu/kAI0eOLPYbjfKkRF9ERCpWSgqWug6Aw9RnfosbuPTSIMd16ADnnQdAHXK4jrd5++2q13daRORUoURfREQq1owZ/sUPuYKhoxtSu6ipIMaM8S/eyJukpcG3mgdNROSkKNEXEZGK4xwu4FH1MxjN1VcXc/w110DdugD0YAVt+J6PPqrgGEVEqikl+iIiUnGysvhvn6tZTRz7acyK5kPye+cEFxkJF1zgXx3KXObNq/gwRUSqIyX6IiJScSIieK7NE8SzmnakMeTKcMLCSihzxRUA7KEZ9TjC8uXw448llBE5hX3wwQeYGRs3bixTPWPHjj3ugU7BPProo8etn1fsJ++iTZ48mSeffPKkyuZLTExk6NChxR6zf/9+/v73v/vXd+7cWW4DdWsCJfoiIlKhPvvM9/MAkVx+eSkKjBwJX33FNf1/YCq+mXg+/rji4hMJtenTp9OvX78Sk/TyUDDR//rrryu8zbIomOifeeaZlfI+VRdK9EVEpMLs2QMrV/qWw8KgVM+IOv10OO88Lrns2K3/Tz+tkPBEQu7QoUN89dVXvPrqq7z33nv+7fkPVRsxYgSdO3fm+uuvJ/8hp1OmTKFXr17ExsYyfvx4Cj789PPPP+eqq67yr3/22WdcffXV3HfffWRmZpKQkMD1118PQMOGDf3HPfHEE8TFxREfH899990HwCuvvEKvXr2Ij4/nV7/6FRkZGcWez7vvvktsbCzx8fEMGDAAgKysLMaNG0dcXBzdunVj4cKFhcoV/IYgNjaWtLQ07rvvPrZs2UJCQgKTJk0iLS2N2NjYYuudNm0aV199NUOGDKFjx47cc889JfwWqi8l+iIiUmE+//zY9JjnnguNG5e+7C9/eWz5iy80zaZUT7Nnz2bIkCF06tSJJk2asGLFCv++lStXMnXqVNatW8fWrVv56quvAJg4cSLLli0jJSWFzMxM5s6de1ydF154IampqezZsweA119/nXHjxvHYY48RERFBcnIyb7311nFl/vOf/zB79my++eYbVq1a5U+Or776apYtW8aqVavo0qULr776arHnM2XKFD755BNWrVrFnDlzAHj++ecBWLNmDdOnT2fMmDFkZWWV6v157LHHOOecc0hOTuavf/3rcfuKqzc5OZmZM2eyZs0aZs6cyfbt20vVXnWjRF9ERCrGvfcS/4fB/D+mcDZbuOiiEyuekACNGvmW//tf+O678g9R5DiTJ4NZ6V7jxxcuP3788cdMnlxik9OnT2f06NEA/OpXv2L69On+fb1796Z169bUqlWLhIQE0tLSAFi4cCF9+vQhLi6OBQsWsLbAI6TNjBtuuIF//etf7N+/nyVLlnBp0IdXHDN//nzGjRtH/fr1AfjFL34BQEpKCv379ycuLo633nqrUFsFnX/++YwdO5ZXXnmF3NxcABYvXswNN9wAQOfOnWnbtm2ZxyOUVO/gwYNp3Lgx4eHhREdHs23btjK3dyoqaiZjERGRsvn4Y7rsWs0UFrCEvlx00TknVLz2zu95sN1Cmq1ZwF08xaJFzTj77AqKVSQE9u7dy4IFC0hJScHMyMnJoVatWjzxxBMA1KtXz39sWFgYOTk5ZGVlcfvtt5OUlESbNm2YPHly0Lvj48aN44orriA8PJyRI0dSu8iHV/g45zCzQtvHjh3L7NmziY+PZ9q0aSQmJhZbz4svvsg333zDRx99REJCAsnJyYW6FgVTu3Zt8vLy/OulueNfXL3B3ruaSHf0RUSk/O3bh1uzBoAcwkgO70vPnidYx+jR3LVmLGN4k4F8waJF5R+mSCjNmjWLG2+8kW3btpGWlkZqairt27dn8eLFRZbJT4CbNWvGoUOHihyYeuaZZ3LmmWfy8MMPM3bsWP/2OnXqkJ2dXej4iy++mNdee83fB3/fvn0ApKen07JlS7Kzswt19wlmy5Yt9OnThylTptCsWTO2b9/OgAED/GU3btzI999/T1RU1HHl2rVr5++2tGLFCr7zvsJr1KgR6enpQdsqTb01nRJ9EREpf4sXY97dtuX0IK5vw/znYJXehRf6Fy9gIV98UY7xiQQzebJvMEhpXi+/XLj8yy8ff0wJXXemT59+3KBZ8HXfefvtt4ssExkZyW9+8xvi4uK48sor6dWrV5HHXn/99bRp04bo6Gj/tvHjx9O1a1f/YNx8Q4YMYdiwYfTs2ZOEhAT/wNiHHnqIPn36cNFFF9G5c+dizwdg0qRJxMXFERsby4ABA4iPj+f2228nNzeXuLg4Ro0axbRp0467455/3vv27SMhIYEXXniBTp06AdC0aVPOP/98YmNjmTRp0nFlSlNvTWel+TpFCouKinIbNmwIdRhSheTPkCASqMZeF3ffDU89BcATTOLwA0/w5z+fYB2ffgqXXAJAEj3oRRI//uiblOdUV2OviyomNTWVLl26hDoMv/T0dBrlD0wpBxMnTqRbt27cfPPN5VanVKxg16SZLXfOneh3ooDu6IuISEX48kv/4iIG0L//SdTRu7d/MZ5VRJDBN9+UQ2wiNUCPHj1YvXo1v/71r0MdioSQEn0RESlfWVm4/MnzgW9r9eXcc0+inshI8Loc1CGHHixXoi9SSsuXL2fRokXqylLDKdEXEZHytWoV5g3220QH2iQ0JeCZPCcm4BPCuSzl22/LIT4RkRpCib6IiJSvgGz8G/rQp08Z6gqS6AfMwCciIsVQoi8iIuUroH/Nt/QO7Gp/4gIS/b4s4cABRzk8Z0dEpEZQoi8iIuXKJSf7l7+lN8XM/ley6Gj/43HPZBdt2E5SUhkDFBGpIZToi4hIudrx7xX0IInb+DubGyRQiqm3ixYWdtzsO+eylIBxviKntKysLHr37k18fDwxMTE88sgjJZaZPXs269atq4ToSufRRx8NdQgnbNeuXQwdOhTwPZ34ggsuoGHDhkycOPG445YvX05cXBwdOnTgd7/7XYlP+B09ejSbNm2qsLhPhhJ9EREpV8tW1WUFPXiR24jrFU5YWBkrHD2aTVf8gZG8w0IuIOALA5FTWr169ViwYAGrVq0iOTmZ+fPns3Tp0mLLVHain5ubW+z+E030nXPkncRAm5ycnGLXS1sO4Omnn+Y3v/kNAOHh4Tz00EP+B4QFuu2223j55ZfZtGkTmzZt4uOPPy62rdtuu40nnniiVHFVFiX6IiJSrgJnxilT//x8t9xCnWefYhYj+YnmrFzpe+ioyKnOzGjoTUmVnZ1NTk4OZgZAu3btuPfee+nduze9e/dm8+bNfP3118yZM4dJkyaRkJDAli1bjqtv27ZtDB48mK5duzJ48GC+//57AMaOHcuECRPo378/nTp1Yu7cuYAviZ80aRK9evWia9euvPTSS4DvgW4XXHAB1113HXFxcQBceeWV9OjRg5iYGF72ngp83333kZmZSUJCgv9Ju08//TSxsbHExsYydepUANLS0ujSpQu333473bt3Z/v27cfFvXz5cgYOHEiPHj245JJL2LVrFwCDBg3ij3/8IwMHDuRvf/sbY8eO5Q9/+AMXXHAB9957L/v27ePKK6+ka9eunHvuuaxevRqAyZMnM378eC6++GJuvPHGQu/7e++9x5AhQwBo0KAB/fr1Izw8/Lhjdu3axcGDB+nbty9mxo033sjs2bPJycmhV69eJCYmAnD//ffzpz/9CYD+/fszf/78Un8IqRTOuZC9gCHABmAzcF+Q/fWAmd7+b4B2Afvu97ZvAC4J2P4asBtIKVDXX4H1wGrgAyDS294OyASSvdeLpYm9U6dOTiTQwoULQx2CVEE18bq44ALnfKm4c+++Wz515uU5Fxl5rN60tPKpN1Rq4nVRFa1bty7UIbicnBwXHx/vGjRo4O68807/9rZt27qHH37YOefcG2+84S6//HLnnHNjxoxx7xbxD2vo0KFu2rRpzjnnXn31VTd8+HB/mUsuucTl5ua6jRs3ulatWrnMzEz30ksvuYceesg551xWVpbr0aOH27p1q1u4cKGrX7++27p1q7/uvXv3Ouecy8jIcDExMe6nn35yzjnXoEED/zFJSUkuNjbWHTp0yKWnp7vo6Gi3YsUK99133zkzc0uWLCkU89F23D0eAAAgAElEQVSjR13fvn3d7t27nXPOzZgxw40bN84559zAgQPdbbfd5j92zJgx7vLLL3c5OTnOOecmTpzoJk+e7Jxz7vPPP3fx8fHOOecefPBB1717d5eRkVGova1bt7ru3bsX2v7666+7O+64w7++bNkyN3jwYP/6okWL/L+DlJQU17lzZ/fpp5+6hIQEd+TIEf9xv/zlL11SUlKh+ksr2DUJJLmTzLVrV/QHiaKYWRjwPHARsANYZmZznHOB30fdDPzsnOtgZqOBx4FRZhYNjAZigDOB+WbWyTmXC0wDngPeLNDkZ8D9zrkcM3sc3weFe719W5xzCRVyoiIiNUjefz7hv9/G4vvTbOVzRx8wg4QE8G6ikZwMbduWT90iAPZnq7C63YNFfwUVFhZGcnIy+/fvZ9iwYaSkpBAbGwvAtdde6//5P//zPyW2s2TJEt5//30AbrjhBu655x7/vmuuuYZatWrRsWNHzj77bNavX8+nn37K6tWrmTVrFgAHDhxg06ZN1K1bl969e9O+fXt/+WeffZYPPvgAgO3bt7Np0yaaNm16XPuLFy/mqquuokGDBgBcffXVfPnllwwbNoy2bdtybpAn523YsIGUlBQuuugiwPctQ8uWLf37R40addzxI0eOJMzrD7h48WLee+89AC688EL27t3LgQMHABg2bBgRERGF2tu1axfNmzcv8b10Qb42zP+2JSYmhhtuuIErrriCJUuWULduXf8xp59+Ojt37qRHjx4ltlEZQpboA72Bzc65rQBmNgMYDgQm+sOByd7yLOA5873Lw4EZzrkjwHdmttmrb4lzbpGZtSvYmHPu04DVpcCIcj0bEZGabvdual02hA1AGm3p0/w72rQpv+SpWzdYlJhLOFmsXNmA4cPLrWqRkIuMjKRfv358/PHH/kQ/P7EsuFxaxZU3M5xz/N///R+XXHLJcfsSExP9yXr++vz581myZAn169dn0KBBZGVlFWovWHKcL7C+gmViYmJYsmRJqcoFrheXjBfVXkRERNDYC2rdujU7duzwr+/YsYMzzzzTv75mzRoiIyP58ccfjyuXlZUV9ANGqISyj34rILCT1g5vW9BjnHM5wAGgaSnLFucm4D8B6+3NbKWZfWFm/U+gHhERyRcwSnYXLenR0ziJ3CS42bO578PzOUBj/sQjGpAr1cKePXvYv38/AJmZmSQmJtI5YJqqmTNn+n/27dsXgEaNGpGenh60vvPOO48ZM2YA8NZbb9GvXz//vnfffZe8vDy2bNnC1q1biYqK4pJLLuGFF14g23uS9caNGzl8+HCheg8cOECTJk2oX78+69evP27AcJ06dfzlBwwYwOzZs8nIyODw4cN88MEH9O9ffFoVFRXFnj17/Il+dnY2a9euLbZMvgEDBvDWW28Bvg8jzZo147TTTiu2TKdOnUhLSyux7pYtW9KoUSOWLl2Kc44333yT4d7dhffff5+9e/eyaNEifve73/l/h+B7D2NiYkoVf2UI5R39YH/+C340K+qY0pQN3qjZn4Ac4C1v0y7gLOfcXjPrAcw2sxjn3MEgZccD4wGaN2/uH4ghAnDo0CFdE1JITbou2rz3Hud4y8kk0KTJNhITvyuXuk9fsYLozV8D0I2VPLs0i8TE4mcnqcpq0nVRlTVu3NifNB/8Q6H/9stNUYn55s2bmTBhArm5ueTl5XHllVcycOBA0tPTcc5x8OBBevbsSV5eHq+99hrp6ekMGzaM3/72t0ydOpU333yTs88+21/fo48+yh133MHjjz9Os2bN+Pvf/056ejrZ2dm0b9+efv36sXv3bp5++mmys7MZNWoUGzduJCEhAecczZo14+233yYjI4OcnBx/3Oeffz7PPfccsbGxdOzYkV69epGRkUF6ejpjx44lNjaW+Ph4Xn31Va699lp69uwJwI033kiHDh3Ytm0beXl5Rb4Pb7zxBnfffTcHDx4kJyeH22+/nbPOOovc3FwOHz7sL5ednU1mZqZ//a677uL2228nNjaWiIgI//keOXKEOnXqFNleu3btSE5O5pxzfH+xYmNjOXjwINnZ2XzwwQfMnj2bzp078+STT3LTTTeRmZnJRRddRL9+/UhLS+Oee+7hww8/pGXLltxyyy3cfvvtvPTSS+zevZu6devSsGHDItsuSVZWVvn+bTjZzv1lfQF9gU8C1u/H14c+8JhPgL7ecm3gJ3xJ/nHHBh7njg2wTQnS5hhgCVC/mLgSgZ4lxa/BuFKQBtdJMDXqurj2Wv9o2fG86GbMKMe616/3172LMxw4540FPCXVqOuiCqsKg3EDHTx40L/ctm1bt2fPnnKpt7gBvDXR+++/7/70pz+Ve71PP/20+8c//lGmOsp7MG4ou+4sAzqaWXszq4tvcO2cAsfMwZecg69P/QLvhOcAo82snpm1BzoC31IMMxuCb/DtMOdcRsD25t7AYMzsbK+urWU+OxGRmiagP80q4kkozykOOnYEr89tC36kBbvUfUdETspVV11Fu3btyr3eyMhIxowZU/KBlShkXXecb/abifjuxocBrznn1prZFHyfXOYArwL/9Abb7sP3YQDvuHfwDdzNAe5wvhl3MLPpwCCgmZntAB50zr2KbyaeesBn3kCNpc65CcAAYIqZ5QC5wATn3L7KeRdERKqJjAzchg0YkIexJSKODh3Ksf5atSA+Hr4+1n0nJaUlgweXYxsiVUhp+pGX1rRp08qtrurilltuKfc6x40bV+51llUo++jjnJsHzCuw7YGA5SxgZBFlHwEKPSvaOXdtEccH/S/HOfce8F7poxYRkUJSUjDvaZcb6cQ5XRuU/Ym4BXXr5k/0u7OCtWsvK+cGRESqFz0ZV0REyi6gH00yCeXbbSdft27+xXhWkZJSAW1IjeOKmRJSpDJVxLWoRF9ERMquMhL9rl39i7GksHatb3SuyMkKDw9n7969SvYl5Jxz7N27l/Dw8HKtN6Rdd0REpJookOgPj6+ANqKjcWaYc3RkE0cOZrFjRzht2lRAW1Ij5D8Uac+ePaEOBfBNrVjeiZ6cOsLDw2ndunW51qlEX0REyiynSxwpSzLozHrW0JW4uApopEED7OyzYcsWapNLZ9azdm2CEn05aXXq1KF9+/ahDsMvMTGRbgFd1ETKSl13RESkzJJve4luJNOAwzTo0JKGDSuoIe8TRAYRtGaH+umLiBRDib6IiJRZfs+dPMJI6Bbs4eXl5NFH+dcDG2lEOh8xVIm+iEgxlOiLiEiZrVp1bLlCBuLm69KFNhd2JA/f3J1r11ZgWyIipzgl+iIiUmaBT6mNr4iBuAFiYo4tr1sH3vT9IiJSgBJ9EREpE/fgZIYkPcRVvE9D0gNnwawQzZrBGWf4ljMyoBwfICoiUq1o1h0RETl5zuH+9ix/yvoZgOiG39O6daOKbTMnh8vab+LQj2v4lt6kpLTj7LMrtkkRkVOR7uiLiMjJ272bWgd8SX46DYmMbY1V4FhcAMaN47Wl0bzDKC7mU9avr+D2REROUUr0RUTk5K1bd2yRaGLjKjrLB6Kj/YtxrFGiLyJSBCX6IiJy8gok+oEDZStMwNO4YklRoi8iUgQl+iIicvJCkejHxvoX41jD+lSHc5XQrojIKUaJvoiInDQXikS/bVuc9+jdZuyl7v4f2b27EtoVETnFKNEXEZGTlpdyLNHf2TiaFi0qoVEzrOBdfXXfEREpRIm+iIicnJ9+Iuwn3630DCI4La5txc+4ky/gq4No1inRFxEJQom+iIicnNTUY4t0ITq2Ev9LCZh5pwupSvRFRIJQoi8iIicnFP3z83Xp4l/UHX0RkeD0ZFwRETk5AwbwdJtnaLA9lS/pzy2xJRcpN7qjLyJSIiX6IiJyUvKiuvD/9nYhw1t/pjLv6LdpgzvrLBZ/35Z1dGFn2lEyMupSv34lxiAiUsUp0RcRkZOSlgYZXpbfvLnvVWlq1cK2beOWKNi40bdp40ZISKjEGEREqjj10RcRkZOydu2x5Urtnx+gc+djy+q+IyJyPCX6IiJyUlJSji3HVmb//ABK9EVEihbSRN/MhpjZBjPbbGb3Bdlfz8xmevu/MbN2Afvu97ZvMLNLAra/Zma7zSylQF2/MLPPzGyT97OJt93M7FmvrtVm1r3izlhEpJp4/XVufvRsPmQo1/K27uiLiFRBIUv0zSwMeB64FIgGrjWz6AKH3Qz87JzrADwDPO6VjQZGAzHAEODvXn0A07xtBd0HfO6c6wh87q3jtd/Re40HXiiP8xMRqdZSUjj90HcM5SPOYUtoEv0jR+i/+XWeYBJvcKMSfRGRAkJ5R783sNk5t9U5dxSYAQwvcMxw4A1veRYw2MzM2z7DOXfEOfcdsNmrD+fcImBfkPYC63oDuDJg+5vOZykQaWYty+UMRUSqKbc2hHPo5wsL45y/3soknuRG/snO9QfJzQ1BHCIiVVQoZ91pBWwPWN8B9CnqGOdcjpkdAJp625cWKNuqhPbOcM7t8uraZWanFxNHK2BXwQrMbDy+u/40b96cxMTEEpqUmuTQoUO6JqSQ6npd9ExKpqG3vLNxB1avTgxJHL1ataJBWhoA7Y6s55138mjZMisksZyI6npdSNnoupDyFspE34Jsc6U8pjRlyzMO30bnXgZeBoiKinKDBg06ySalOkpMTETXhBRULa+L9HTY+wMA2dSmcY9oBg2qE5pYevb0zfOJ7wm5kZFjORXe7mp5XUiZ6bqQ8hbKrjs7gDYB662BnUUdY2a1gcb4uuWUpmxBP+Z3yfF+7j6BOEREJF9AZ/iNdKJzXIiSfNATckVEihHKRH8Z0NHM2ptZXXyDa+cUOGYOMMZbHgEscM45b/tob1ae9vgG0n5bQnuBdY0B/h2w/UZv9p1zgQP5XXxERCSIdcf3zw/V1JoAdOniX4xmHampIYxFRKSKCVmi75zLASYCnwCpwDvOubVmNsXMhnmHvQo0NbPNwB/wZspxzq0F3gHWAR8DdzjncgHMbDqwBIgysx1mdrNX12PARWa2CbjIWweYB2zFN6D3FeD2CjxtEZFT37oqMBA3X4E7+hs2hDAWEZEqJpR99HHOzcOXaAdueyBgOQsYWUTZR4BHgmy/tojj9wKDg2x3wB0nFLiISA2Wl7LOf5doHdHcWXBi5MrUqROuVi0sL4+z2cq29ZlARAgDEhGpOvRkXBEROSE5a47d0d93RjSNG4cwmPBwaN8egFo4Gu/eyP79IYxHRKQKUaIvIiKll5lJnR3fAZBLLSLiO4U4IDB13xERCUqJvoiIlF5EBI/d8zPnsoTreJtOcfVCHVGhAblK9EVEfELaR19ERE49K7Y05hvO5RvO5bVQDsTNd/HFfLkoj9eWRvM15/ErTbEpIgIo0RcRkROUknJsOaQz7uQbPJi0nYOZ5j0vXXf0RUR81HVHRERK7cgR2LTp2Hp0KGfcCdC587FlPTRLRMRHib6IiJTa9lnf0CR3DwDt2kHDhqGNJ19U1LHlzZshNzd0sYiIVBVK9EVEpHSOHOHsG89nD6fzA2cQ3+VoqCPyO+00aNECwJF7NIe0tBAHJCJSBSjRFxGR0tm0iVp5vlvlh2hIVFzdEAcU4J13+CSzP3tozj08oe47IiIo0RcRkdJad+xBWeuIJjY2hLEUtH8/XQ8sphl7NcWmiIhHib6IiJROgUS/Ssy4k09z6YuIFKJEX0RESiVnzbFEP5Xo42a6CbmA6X86s54NqXkhDEZEpGpQoi8iIqWSvepYon+gVTT164cwmIKaNiWn6ekA1CeTjNRtIQ5IRCT0lOiLiEjJsrOpm7bRv1ovvirdzvcJiz12V7/5T+vYvz+EwYiIVAFK9EVEpGRbthCWmw3ANs7inIRGIQ6oMIs5luirn76IiBJ9EREpjao8EDefBuSKiBxHib6IiJSsKk+tmS9ad/RFRAIp0RcRkRIdqdOAVDqTQxjrLZqoqFBHFERAot+FVNanuhAGIyISekr0RUSkRGsG/w/RpFKfDL7tcB316oU6oiDOOIOc05oAEEEmP639McQBiYiElhJ9EREpUUqK72c2dekQFxHaYIpiRtY/ZxHNWhpwmKVpLcjNDXVQIiKho0RfRERKtHbtseUq2T/f03DYhfzcIpps6nL0KKSlhToiEZHQUaIvIiIlCkz0q+SMOwECxw9oQK6I1GRK9EVEpHgffUS/JU9wOXM5nR+rfKLfOeBZXkr0RaQmC2mib2ZDzGyDmW02s/uC7K9nZjO9/d+YWbuAffd72zeY2SUl1WlmX5pZsvfaaWazve2DzOxAwL4HKvasRUROLUffnM4f99/LXK5gRK0P6Ngx1BEVr3PHXM5mC5fxkWbeEZEarXaoGjazMOB54CJgB7DMzOY459YFHHYz8LNzroOZjQYeB0aZWTQwGogBzgTmm1knr0zQOp1z/QPafg/4d0A7XzrnhlbMmYqInNqyV66lrrecflYMdesWe3hoOcdtD5/JnewG4OqUXUCL0MYkIhIiobyj3xvY7Jzb6pw7CswAhhc4Zjjwhrc8CxhsZuZtn+GcO+Kc+w7Y7NVXYp1m1gi4EJhdQeclIlJ95OZSL229f7V2fBXvt2OGa9fev1pr/bpiDhYRqd5Cmei3ArYHrO/wtgU9xjmXAxwAmhZTtjR1XgV87pw7GLCtr5mtMrP/mFkV/19MRKQSbd1K7ewsAHbRgnbdfxHigEpWN76Lf7nFz+s4cCCEwYiIhFDIuu4AFmRbwc6URR1T1PZgH1wK1nkt8I+A9RVAW+fcITO7DN+d/qA9UM1sPDAeoHnz5iQmJgY7TGqoQ4cO6ZqQQk7166LZ4sXkz6a5lhgghcTEn0IZUonahNfjHG+5C6m8/fZyunRJD2lMBZ3q14VUDF0XUt5CmejvANoErLcGdhZxzA4zqw00BvaVULbIOs2sKb7uPVflbwu8s++cm2dmfzezZs65Qv+TOedeBl4GiIqKcoMGDSrViUrNkJiYiK4JKeiUvy6++sq/mEIso0bFHjd9ZZV0+DC89BIA0axje4MeVLVfwSl/XUiF0HUh5a3ErjtmVt/M/p+ZveKtdzSz8hi4ugzoaGbtzawuvsG1cwocMwcY4y2PABY455y3fbQ3K097fHfgvy1FnSOBuc65rIDza+H1+8fMeuN7T/aWw/mJiJzyjq48NoH+hrAYzjmnmIOriujoY4us0xSbIlJjleaO/uvAcqCvt74DeBeYW5aGnXM5ZjYR+AQIA15zzq01sylAknNuDvAq8E8z24zvTv5or+xaM3sHWAfkAHc453IBgtUZ0Oxo4LECoYwAbjOzHCATGO19mBARqfGyk4/NuHO4XQy1Q/k9cGm1bUtO3QhqH83kDHazc/VPQLNQRyUiUulK8yf7HOfcKDO7FsA5l5l/B7ysnHPzgHkFtj0QsJyF7y58sLKPAI+Ups6AfYOCbHsOeO5E4hYRqRFycgjfdmzGnYge0cUcXIXUqsWR9p2pvWElALkpqUD/4suIiFRDpZl156iZReANajWzc4AjFRqViIiE3pYthOUcBWAHrejQMzLEAZVenfhjH0pO27GO3NwQBiMiEiKluaP/IPAx0MbM3gLOB8ZWZFAiIlIFNG7Mc2c9QcPv13KQ0+jaNdQBlV7d+Gh4x7fcISeVbdvg7LNDG5OISGUrMdF3zn1mZiuAc/FNa/n7YDPSiIhI9ZJ3egvu3zeJQ976rviQhnNioqM5HNaIlNwubKcN69cr0ReRmqfIRN/MuhfYtMv7eZaZneWcW1FxYYmISKilpcEhL8tv3hzOOCOk4ZyYK67gDzcd4OVXfEPKWm+Ayy4LcUwiIpWsuDv6T3k/w4GewCp8d/S7At8A/So2NBERCaVVq44td+0K5TMNQyUJC6PzsQfkaopNEamRihyM65y7wDl3AbAN6O6c6+mc6wF0AzZXVoAiIhIaq1cfWz6V+ufnC3yw1/r1RR8nIlJdlWYwbmfn3Jr8FedcipklVGBMIiISauvXc8uTw4kjhkUMoGvXO0Md0QkLTPR1R19EaqLSJPqpZvYP4F/4ptj8NZBaoVGJiEhorVlDq0MbuZqNhJNFy/hTL9FvV383vw77nC65a9j+QxsOHLiNxo1DHZWISOUpTaI/DrgN+L23vgh4ocIiEhGRkDu68tgTcVMthgu7FHt4lRSWvJx/5l4HwGLOZ8OG2+jdO8RBiYhUotJMr5kFPOO9RESkBkj/Zi1NveV9LWMIDw9pOCcnLs6/GEsKc1IdvXufSiOKRUTKpsQn45rZd2a2teCrMoITEZHQqLUu5dhKTEzoAimLVq3ICG8CQCQH+O/S7SEOSESkcpWm607PgOVwYCTwi4oJR0REQi4zk8Y/bgQgD+MX/aJDHNBJMuNQ+zjqpy4C4EjSGuCs0MYkIlKJSryj75zbG/D6r3NuKnBhJcQmIiKhsHYttVweAJvpQJeeDUIc0MkLiz/Wfaf+5tXFHCkiUv2UeEe/wBNya+G7w9+owiISEZGQcsmryO/Jvop4zo0r9vAqrXH/rjDDt9x6/xoOHYKGDUMbk4hIZSlN152nApZzgO+AayomHBERCbUDi1YR6S1vrh/PiNYhDadManc79ikljjWsW4dm3hGRGqM0if7NzrnjBt+aWfsKikdERELsaNIq/3JWVDx2Kk9UExvrX+zMet5KPkrv3nWLKSAiUn2U2EcfmFXKbSIicqpzjkZbj/Vlr983PoTBlINGjfg5sh0Adcjhp8XrQxuPiEglKvKOvpl1BmKAxmZ2dcCu0/DNviMiItWNGdf33crexNXEksIFg9qEOqIyy+zQlSZJaQDkJq8BuoY0HhGRylJc150oYCgQCVwRsD0d+E1FBiUiIqHhHHyxugn7GMgiBnJXj1BHVHZ25XD+L+ksVtOVNbvP595QByQiUkmKTPSdc/8G/m1mfZ1zSyoxJhERCZHvv4d9+3zLkZHQvhqMyGp+703cPQWOHgV+hP37fecmIlLdFdd15x7n3BPAdWZ2bcH9zrnfVWhkIiJS6VasOLbcvTun9kBcT+3a0LkzrPaGHqxbB+edF9qYREQqQ3GDcVO9n0nA8iAvERGpZn6cs5TG7Ad8iX51ERNzbDklJXRxiIhUpuK67nzo/Xyj8sIREZGQOXSI8dPOYwKOLZzNsvgNlG4W5qrv2CybjvWrswFNsSki1V9xXXc+BFxR+51zwyokIhERCQm3eg21vD/7WYTTvXf1SPIBLjz8IR/zHN1ZwUfzfgs8EOqQREQqXHF/xZ+s6MbNbAjwNyAM+Idz7rEC++sBbwI9gL3AKOdcmrfvfuBmIBf4nXPuk+LqNLNpwEDggFf9WOdcspmZd/xlQIa3PaCXqohIzRD4RNy1teMZ0SGk4ZSrsyP3cTqfAnDGTv2JF5GaobiuO1/kL5tZXaAzvjv8G5xzR8vasJmFAc8DFwE7gGVmNsc5ty7gsJuBn51zHcxsNPA4MMrMooHR+Ob5PxOYb2advDLF1TnJOVfwYV+XAh29Vx/gBe+niEiNsv+LY4n+vtbx1CrNIxVPEc0u7g73+JZjjqxgzx5o3jy0MYmIVLQS/4yb2eXAFuBZ4Dlgs5ldWg5t9wY2O+e2eh8cZgDDCxwzHMgfIzALGOzdgR8OzHDOHXHOfQds9uorTZ0FDQfedD5LgUgza1kO5ycickqpszrJv5yXUI1G4gK1YrqQZb5nPZ7FdtZ9sSfEEYmIVLzSdMB8CrjAObcZwMzOAT4C/lPGtlsB2wPWd1D4Trr/GOdcjpkdAJp625cWKNvKWy6uzkfM7AHgc+A+59yRIuJoBewqGLCZjQfGAzRv3pzExMQST1JqjkOHDumakEJOlevCsrPpu2u1f31vu4anRNwnomXjLkTtXwnA6jfm4pqF7iEBp8p1IZVL14WUt9Ik+rvzk3zPVmB3ObQdbHbmgoN/izqmqO3BvqHIr/N+4Ad8Uy28DNwLTCllHL6Nzr3slSUqKsoNGjQo2GFSQyUmJqJrQgo6Va4Lt3wF5vXK3Ep7Rt3Wl06dSigUAlt/3srjix9n075N9G7Vm3vPv5cmEU1KVTYltg8s9iX6p+/4kUGDxlVkqMU6Va4LqVy6LqS8lSbRX2tm84B38CXAI/H1fb8awDn3/km2vQNoE7DeGthZxDE7zKw20BjYV0LZoNudc/l36I+Y2evA3ScQh4hItbbnP0mc7i2vrtOT4R1DGk5QyT8kM3DaQA4eOQjAwrSFvLvuXRaPW0zLRiX3uAw/rzss9i033rqyIkMVEakSSjPUKhz4Ed+MNYOAPcAvgCuAoWVoexnQ0czae4N9RwNzChwzBxjjLY8AFjjnnLd9tJnVM7P2+AbSfltcnfn97r0+/lcCKQFt3Gg+5wIHAj4UiIjUCAcWHHsO4t72PavcE3EzsjO4auZV/iQ/39aft3Ld+9fh+6+heC0vPzbu4JyDK8jKKvcwRUSqlBLv6DvnKuS7Ta/P/UTgE3xTYb7mnFtrZlOAJOfcHOBV4J9mthnfnfzRXtm1ZvYOsA7IAe5wzuUCBKvTa/Its//f3p3HR1Xd/x9/fSYbISSEJezIIpvILqLWBXBFq6XuW1FbqlZrW2tta/X3bfu12qpVW7uItdWvu7iCKFIrloArAoKy76BAkD2EJGQ9vz9mMpnJMpkksyW8nz7mwbl3zj33M3Cd+cyZc8+xHLxDdZYBP/Dtfxvv1Job8E6vGb/fckVE4qTNiuobcVNOPC6OkdTtoY8eYsuBLQBkpWVx89ibeeCjB6h0leRuyWX6iulcOfzKkG1knDCMMpJJoZyBbGDZwnxGjW8fg+hFROKjwUTf12P+I6BvYP1ILJjlnHsbb6IduO/XAeXDeIcK1XXsvcC94bTp2396Pe044IeNClxEpDVxjuWVx1JCPv3ZRI/zE2vGncLSQv688M/+7QfPepDrj7ue8spyHvzYu+TL3Qvu5vJhl+OxED9Up6Wxrf0w+uUvA/CaIVoAACAASURBVGDHnGWMGj8+qrGLiMRTOEN3ZgJbgL/inYGn6iEiIq3A4RLj2wefYSAb6MB+xpwR3s2tsfLC8hfYV7wPgH7Z/fjuaO8Pr3eddheZqZkArNmzhnc2vNNgW/lHe7/ElJDK3qVfNVBbRKRlCyfRP+yc+4tzbp5zbn7VI+qRiYhITCxbBmVl3nK3Qe3p2DG+8dT09OdP+8s/Gvcjkj3eH5ez22Rz/Zjr/c8988UzDbZ14Ls/ZRwLyeIgT5R8J/LBiogkkHAS/UfM7DdmdpKZjal6RD0yERGJiU8CViU5IcHWBd9+cDsffvUhAMmeZK4ecXXQ89eOutZfnrlmZq2bdWs6evIwFjGOUtJYtgzCuIdXRKTFCifRHw5cD9xH9bCdB6MZlIiIxM5HH1WXEy3Rn71+tr88oe8EumR0CXp+RNcRjOg6AoDD5Yd5fXXoGZ979cL/i0V+PmzeHNl4RUQSSTiJ/oVAf+fceOfcRN+jzhtbRUSkZXH79nPNrEu4jYcYx0JOOSXeEQV7a91b/vL5A+ue0fk7w6uH4Ly88uWQ7ZnB2LHV24sWNS8+EZFEFk6i/zmQHe1AREQk9vJmfML5Ja/xELczLekWhg2Ld0TVisuKmbtprn/7m4O+WWe9S4Ze4i/P2zKP4rLikO2OHQvdyGMyMyl4eU5kghURSUDhJPpdgTVm9o6ZzfI93oh2YCIiEn273vjYX97W6ySSkuIYTA25W3IpLvcm7YM7DWZAxwF11uvXoR9DOg8BvMN3crfkhmz3fJtNHj2YyYWMfe++iMYsIpJIwkn0f4N3+M7vgYfxrkBb97utiIi0KClLqhN9d+JJcYyktvlbqyd4O3fAuSHrBj4/Z0PoXvq+F472lwflL6KipLyJEYqIJLYGE33fVJr5wDeBp4AzgMeiG5aIiERdRQV98hb6N3tekliJ/vtfvu8vn9bntJB1zxt4nr/89vpaayYG6TamB9s9vQFoSzFb31rejChFRBJXvYm+mQ0ys1+b2Wrgb8BXgPluxv1rzCIUEZGo2LNgFe1cAQB5dGP4+X3iHFG14rJiFm2vvlP25KNODln/1KNOpW1KWwA27t/I1gNb661rBpu7nujf3jXzo3rrioi0ZKF69Nfg7b2/wDl3ii+5r4hNWCIiEm1bplcP21nf6STS2lgcowm2eMdiyiq9q3gN6jSo1rSaNaUlp/GN3t/wby/YuiBk/cJR1V8cUha+H6KmiEjLFSrRvxjYCcwzs3+a2RlA4nwKiIhIs1TM/9BfLhqZWMN2PvjyA3/51KNODeuY8X3G+8uB4/vr0u686qFAR21ZoJWzRKRVqjfRd87NcM5dDgwBcoGfAl3NbJqZnR2j+EREJEp6baxOhjtODi+ZjpWq1XABTjkqvMn9G5PoD750BAdoD0BOWR6lqzc2IUoRkcQWzs24hc65551z5wO9gGXAHVGPTEREoibv4y30LPeOYz9EBsOvOy6+AQVwzrF4x2L/9om9TgxRu9q4nuNok9wGgA37NrCjYEe9dTt3TeKz9OovENteCP3FQESkJQpnek0/59w+59w/tDKuiEjLtvGJXH95dcdTSM9KiV8wNewo2MHXhV8D0C61HYM6DQrruLTktKAvBe9vDT32fteQ6uE7xe+EHtMvItISNSrRFxGR1uGJwis4jfn8mv9lx9nXxTucIIG9+aO7jcZj4X9Undy7+ibbT7Z9ErJu2lnViX7nVUr0RaT1SY53ACIiElvOwZx5bfia03if01j0s3hHFGxJ3hJ/+bjujRtSFNijv3D7whA1YeAVx/H5AyNYzFiWJo3nrxWVWJL6v0Sk9VCiLyJyhFm5Er72joyhY0cYPTp0/VgLSvR7NC7RP6HnCf7yZ3mfUVpRSmpSap11h45MoVP25xw4ABTATzbBwIFNCllEJCGp60JE5Agzd251+fTTISkpfrHU5JxjyY7qRH9sj7GNOj4nI4d+2f0AKKko4Yuvv6i3rscD36ieep8PPqi3qohIi6REX0TkCLP5hY/pxB4AzjwzzsHU8HXh1/4bcTNSMsK+ETdQ4PCdhsbpnxIwc+eHH9ZfT0SkJVKiLyJyBDl4oJI7F32bXXThE07g/LE74x1SkOVfL/eXh3cd3qgbcasEDt9paJz+ydX37vLJ+2VQWNjo84mIJCol+iIiR5BP/7GUruzCg2NQ0iZ6ju4S75CCrNi1wl8eljOsSW0E3ZC7LXSif/zxcE7SXGbwbT5a14mC+x9t0jlFRBKREn0RkSNI/vQ5/vKXx5zjHaieQJbvCu7Rb4pR3Ub5b8Bdv289e4v21ls3PR3G9/+Kb/MGWRRQ+Po7TTqniEgiius7vJlNMrO1ZrbBzGqttmtmaWb2ku/5hWbWN+C5X/n2rzWzcxpq08ye9+1fYWZPmlmKb/8EM8s3s2W+x6+j+6pFROKjogJ6rahO9DMvOzeO0dQtqEe/S9N69NOS0xjVbZR/+9Ptn4as3+aCs/zlTqvfh6KiJp1XRCTRxC3RN7Mk4O/AucBQ4EozG1qj2lRgv3NuAPAn4H7fsUOBK4BjgUnAo2aW1ECbzwNDgOFAOvD9gPO875wb5XvcHflXKyISf4ve2cfYcu/NqZUYfW88p4EjYqvSVbJy90r/9vAuTevRh+Bx+ot2LApd9+JerOBYAFIqS3G585t8XhGRRBLPHv1xwAbn3CbnXCkwHZhco85k4Glf+VXgDDMz3/7pzrkS59xmYIOvvXrbdM697XyAT4FeUX59IiIJZcNjc0miEoCtOcfj6dI5zhEF27x/M0Vl3t70LhldyMnIaXJbx/c43l9uKNEfNw7mpVZ/6cl/dlaTzysikkjimej3BL4K2N7m21dnHedcOZAPdApxbINt+obsTAH+HbD7JDP73MzmmNmxTX1BIiKJyjlI/+9b/u2yMybFMZq6RWLYTpXA+fcX71iMt4+nbsnJkHd8dT9TytszobKyWecXEUkE8VwZ1+rYV/OduL469e2v64tLzTYfBRY45973bX8G9HHOHTKz84CZQJ1rI5rZDcANADk5OeTm5tZVTY5Qhw4d0jUhtSTKdbH683SuLKzuqd5/So+EiCvQrK3V8XUo7dCs+CpdJW2T2lJUUcTOQzt59T+vkpNW/y8ERaO7sevDHLqwm4yDO/ns0Uc5OKx5XzZCSZTrQhKLrguJtHgm+tuA3gHbvYAd9dTZZmbJQHtgXwPH1tummf0GyAFurNrnnDsYUH7bzB41s87OuT01A3bOPQ48DjB48GA3YcKEsF6oHBlyc3PRNSE1Jcp1seZPc8gmH4A97fpwws03gNXVZxI/j736mL88acwkJoyZ0Kz2xn05jtwtuQAkH5XMhGPqb697d3jjb5O5nn8BMHLzVpJuuaVZ5w8lUa4LSSy6LiTS4jl0ZxEw0Mz6mVkq3ptraw6MnAVc6ytfAvzXN8Z+FnCFb1aefnh74D8N1aaZfR84B7jSOef/TdbMuvnG/WNm4/D+ndQ/F5uISAtTUQFZ777q3y4895KES/KhxtSazbgRt8rY7sHDd0IZNAg+6Hyhf7tk+gzveCcRkRYsbj36zrlyM7sFeAdIAp50zq00s7uBxc65WcATwLNmtgFvT/4VvmNXmtnLwCqgHPihc64CoK42fad8DNgKfOzL61/3zbBzCXCTmZUDxcAVLtRgThGRFmb+fPhn8XcowrjIZtDrJ5fEO6RaSitKWbd3nX97aE7NSdga7/ie4d+QawYdLjmDg49lkkUBbXdshBUrYHjzv3CIiMRLPIfu4Jx7G3i7xr5fB5QPA5fWc+y9wL3htOnbX+drdc79DfhbowIXEWlBnnsOcplILhNZcdM0/nxSUrxDqmXDvg2UV5YD0Kd9HzLTMpvdZuDMO1U35FqIXzLOvziN2Y99k4t5jffSz2eSeeq8IUxEpKVIrCURRUQkovLz4aWXqrevmJKScKvhAqzZs8ZfHtJ5SETa7Jvdl07pnQDYf3g/m/ZvCll//Hi4P/NeupPHecWvs7RUk7CJSMuWeO/2IiISMc8/X73Q6/DhcMIJoevHy9o9a/3lwZ0GR6RNMwuaZrOh4TspKTB8cn/24f1yMHNmRMIQEYkbJfoiIq2Uc/DvR9bSngMA3HhjQt6DC8DavdWJfqR69CF4Pv1F20Mn+gAXVt+Py4wZEQtDRCQulOiLiLRSn34Kt6+7njy680LSd5hy6pZ4h1SvwKE7gztHpkcfaozTzws98w7AOedAmzbe8ooVsOaLUigpiVg8IiKxpERfRKSVevX36ziN90nnMJdVTiera3q8Q6qTcy5qPfqBM+8s2bGEisqKkPUzMuC882AYy/kLP6L3iT28459ERFogJfoiIq3Q1q3Q5c0n/NsHTz0funaNY0T121W4iwOHvcOL2qW2o3u77hFru0dmD397hWWFQb8c1GfKFDiHd/gRfyOjeC/ukUc0p76ItEhK9EVEWqG//KGQqe6f/u0Ot0+NYzShBfbmD+40OOQUmE0R2Kvf0MJZ4O3Rn5H9PQppC4B98QXk5kY0JhGRWFCiLyLSyuzaBeVPPE1H9gNQ1ONob/aaoAJn3InksJ0qgeP0G5p5ByA1FSZd1ZGn/QuzA3/6U8TjEhGJNiX6IiKtzCMPlfPD8j/7t9N/+RNISrxFsqoE3Ygboak1AwXOvBNOjz7A9dfDI/zEv+3eegvWr494bCIi0aREX0SkFfnqK9jxp+kMwpuUlmZkY9/7bpyjCi1o6E4EZ9ypEpjoL9u5jNKK0gaPGTUKOp00mNl4fwkx5+DeWouxi4gkNCX6IiKtyK/vLOeOst/5t5NvvxXatYtjRA2Lxqq4gTq37Uy/7H4AlFSUsGLXirCOu/lmuJ9f+rfds8/CqlURj09EJFqU6IuItBKffQbJzz3FYNYBUJ7RHs+tP2ngqPgqKS9h84HNABjGwI4Do3KepgzfufRS2ND9NP7NOd74Kivhf/4nKvGJiESDEn0RkRaspLyEWWtnceuc25j45Dm8/71fceHl8LvT4MtffBeys+MdYkgb92+k0lUC0Ce7D+kp0ZnrP+iG3DBWyAVIS4PbboM7+X31ztdfh4ULIx2eiEhUKNEXEWmB9hTt4X/++z90f6g7k6dP5pFP/8TBnP+w9qg9zDwGfn06HO3+zJQZU9hbtDfe4dYr2jfiVgmcYjOcmXeq3HgjbOkwhpe5tHrn9OmRDE1EJGqU6IuItCBlFWXc/8H99P1zX+55/x72H94fsv5zXzzHiMdGhD0uPdYCp9aMZqI/pvsYDO/8/Ct2raC4rDis4zIz4fbb4Q7uYytHcXP75zl098NRi1NEJJKU6IuItBALty1k7D/Hcsd7d1BYVujfn1zYBz74BUx/ncEf5/LkBU8zacAk//M7CnYw4akJfPH1F/EIO6Q1e6N7I26VrLQs/4w+Fa6CZTuXhX3srbdCSY/+DGAD0/Kv4p57I7ugl4hItCjRFxFJcPmH87nl7Vs46YmTgpL1Y3OOZcym6Ux58E6y5t5Jmy0X8tpD4/numGuYc/Uc3rjiDbLSsgDYW7yXydMns694X7xeRp2CevSjMLVmoMAbchszfKdtW7jnHignBYAHH4Rl4X9PEBGJGyX6IiIJyjnHa6teY+ijQ/n7or/jcAC0TWnLA2f+kdPXLWPMMwd50t3IUkYz467FHHts9fHfGvwt5k6Z60/2txzYwpQZU3DOxePl1OKcC5pDP5o9+hB8Q264M+9UufZaOO00b7miAq65Bop35nuz/gT5+xQRqUmJvohIAtp6YCuTp0/mklcuYUfBDv/+cwecy4qbVrLz9dvZ+Mg7TOMmAPqzmUlrH6nVzvE9j+eZbz/j3357/du8sPyF6L+AMOwq3MWBwwcAaJfaju7tukf1fEEz7zSiRx/A44HHH/fOxANQvnwV+UPGwc9/7h3Er2RfRBKQEn0RkQRSVlHGAx8+wNBHh/Lmujf9+7tmdOWlS15ixiWzue+Ovmx6eAYzuJBkKgBwY8bAtGl1tjl5yGR+PO7H/u3b/nMb+4tD38QbC0Er4nYajFl0x76P7DaSJEvynnvPWg6WHGzU8YMHwyO+71Lf40m65XvXK+Dhh+GXv1SyLyIJR4m+iEgCcM4xY/UMhk8bzi/n/pKisiLAu4jUjcfdyJpb1nBqx8s48wywxx/jFS4llTLvsX37Ym++GXIF3HtOv4eemT0Bb0/6Hz/6Y/RfVAMCx+dHe9gOeIc8DesyDACH47O8zxrdxg03wFVXwa/4A69zYfUTf/wjXHcdFIc3m4+ISCwo0RcRiSPnHHPWz+EbT36Di16+KKiXe2TXkXw09SOmffMxXn8hmxOG5HPTh1fzGDdV9+QPHIQtWAA9eoQ8T2ZaJg+e/aB/+5GFj7CrcFd0XlSYYjWHfqCmLJwVyAyeeAJOPCWFK5jOTCZXP/nMM3DKKbBqVSRCFRFpNiX6IiJxcLDkIP/67F8MmzaM8144j0+2feJ/ListiwfPepBF1y/mwIoTOWFsBR9MfZJFBwdxFS/667lRo7AF86F377DOedmxlzGi6wgAisqKeODDByL7ohopaOhOlGfcqRI4887ivMbdkFulTRuYNQtGjk3lMl7mSb5b/eRnn8GoUXDXXXDgQHPDFRFpFiX6IiIxcrDkIDNWz+Cq166i24PduP7N61m1u7r3NzUpldtOvI0PL91E8qKfMXJ4MueeCxs/O8DfuIWuBPTA33AD9tFH0K1b2Of3mIe7J9zt3358yeONHqceSYE9+rEYugM1VshtQo9+lQ4d4L334NTTU5nKE9zM3ykj2ftkWRn8/vfQp4/3G4GISJzENdE3s0lmttbMNpjZHXU8n2ZmL/meX2hmfQOe+5Vv/1ozO6ehNs2sn6+N9b42Uxs6h4hIUznn+DL/S+bvns8dc+9g/FPj6fRAJy56+SJeXPEixeXVY7kzUtoxuctt3FC4hp23XsX3B6zn1lth9Wrv8/voxDNJ3l7jyp694JVX4B//gPT0Rsd1weAL/El1QWkBTy17qtmvtSlKykvYfGAz4L0PYWDHgTE577Auw0hL8k6ds/nAZvYU7WlyW1lZ8M478POfG9O4mTF8xod8o7rCwYPMXdub/PyAg5yDwsJabYmIRENyvE5sZknA34GzgG3AIjOb5ZwLHNw4FdjvnBtgZlcA9wOXm9lQ4ArgWKAHMNfMBvmOqa/N+4E/Oeemm9ljvran1XeO6L56EWlJKiorKC4vprisOOjPQ6WH2F24mz1Fe9hdtJtdhbvZvG8rG/dvZPOBjRSVh07oeuzrzrilgzl/UQfGHp7LYP5OG0p4j9M5k/cA7/21110H3778pzA3B8/Pfw4ZGU1+LR7z8JMTfsJNs73Tcj6y8BFuGXcLHottv8/G/RupdJUAHNX+KNJTGv+lpSlSk1IZ2W0kn27/FIAlO5ZwzoBzGjiqfsnJ8MADcNFFcNNNwzl12ftczfPcye+pIImzfjEKfgFDhsDYsXBi1mpunjaMkR17kT98DPTtQ3KfXqQd3Yvkrp283x6qHh06NOvfWkQkbok+MA7Y4JzbBGBm04HJQGCiPxn4ra/8KvA3886/NhmY7pwrATab2QZfe9TVppmtBk4HrvLVedrX7rT6zuEaWFFm84F19P5p6L8+Bxz2pHPA0zFof2blQTIqC7x1GphNrsgyKPBkBe3LqjxAG1d7Zoe62iq0TAot4IPCHNmV+0l1JfXHHdBOgbXnsLUJer5D5V6SKSecieTyPe0p8/544texcg8WcHSodvI9HaggyReXwxx0dPX3wNVs64CnI5VUvyAPlbR34U0r6IADng5BO5Ipp50rqN5V8+/839XFSjwUWFZQVMmU0dYV1X0+C46+gmSKfP92Vc+kUkIbdzh03L6YykmhmOB/uzRKSKU0qM36lJFCSY1/uzRKSKY8+Hwhji+zlKBaaZTgobL+2GscX/VvH3h+rOErz+F9/YH/9oD/tfvrhfj/r5xkHOZtLam8/oqNNDoPzl0Pl62EkV/nAXm16pzEx0w6vZRvX5bKVVdBZibAADjltxGJYcqIKdz53p3sP7yfTfs3MW/zPM7of0ZE2g5X4LCdY3KOiem5j+9xvD/RX7RjUbMS/SonngiLFsGrr3p4+OEpDFt0FX3YCr5rcM0a7yOb9zAcHfZ+BblfhWzzjZRLuKHDKyQng0sppCJrI+cl38uIpFnsaethX7rjQLqjKMVRkgwlyY7SJF85CfYmd2BPUhff/zMOhyOncidtXQHOvP+fVPr+dOYP1e+gtafQkxm0r2PFbu//h2HY7+nIYWsbtC+nYmet95D67PXkUGppQfu6VmzHE9anD+xK6kZFQJpjOLpVbA/rWICdST197wFeSZTTpWJnWMdWYnyd1DNoX6oroVPl7pDHOeewN4xyktmdFDwsr40rokNleCtbl5DGvqScoH0ZlQVkufx6jghWbG3ryF0O0M4dCuv4Q9aOAk920L7syn2k1/P5V5OuveBrrznimej3BALf5bYBJ9RXxzlXbmb5QCff/k9qHFv1f1RdbXYCDjjnyuuoX985Qv6eW5YE27IrGniJAId8j2oHfI/wHPQ9qjVu9uvaZwvvbaL+s+1t5vGN+6E8+GyumcdX0NjXH1y7jMb8/VfUql0GhPc2C1Beq3ap7xGeMt+jWonvEc/jwxeJ44OF/3cHhPmBENLh9pA3BrYfT+cdfVj+5Q/p1sDnZFGnXiRPPJU5j+ZDTk7oyk2UkZrBlBFT+MunfwHgiaVPxDXRH9IpNuPzqwTOvLNw+8KItZucDFdcAZdfDp9/nsTrr/dn9mz4/HPvaroAg1hHBR6SQnzhLfPAR73hhaPWsqv3+dB9CWR6E8ynGhXRLt+jWu2vlaHspeZ7aOPmaapdu3Hnr107/FQJvB/1wUJ/tappa9BWeaOP31Lr+PDS3KrawcfXziYaOj74F8V8GvP5U3fuEf7nX+3co/Gf3br2IiGeiX5dX1VqflWqr059++v67TlU/XDj8FY0uwG4AYDoLuAoIonEGZSl06m8iDblkF4G6eWQUQqdiyCnCHIKvX92L4Cj98OV5csppj8pyRWkp8NRbXdSXtGfze3TKWuTzuGsjlR27UBSr/Yk986m+KjeFPXpQ3mmrxdr5cqovqThFcP95VdXvsqVWVeSmZIZ4ojIyl2T6y/bPiM3N7feupFmRdVv+/M2zuO9ee/5F9KKpNNP9z4OH/awYUM7Nm1qxydf386cvNtpu3UTPQ9tpWNRHp1LdtKlbDvbe2/lvTFbeH/wPvLbVgLLfQ8RkaaJZ6K/DQicE64XsKOeOtvMLBloj7eLNdSxde3fA2SbWbKvVz+wfn3nqMU59zjwOEDf3n3ci994vsEX6dLScZ06B+3zFORjhdXfy833XaOuRSFd2wwqs7L9dQA8+fuxkroXZam5sqRrlwntgj+8PQf2Q2ntvtHAQ6vO5zLb49KDf/5K2r8HKiuC6tXHtc/G0oKHj3j27q61gmR9K2K67A6QnFJdD7B9ofsFqtoyA9ehk3ft+ioVFVj+gVrnq3dBzo7B/3ZWXoYdKgje5/s7WLJkCccdd1z1/iQPrn12UB3KSrGi4D6doL/3gA3n8UBm8LAtKy3B6vi3q8kwb/di2xrje0sOQ2l1v7bHE+LfLzkF0tODr8+Sw1Beu6c78Drwv4TUVO8jsN7hEqisrF03sE7VztRUSEoOrldyGJwLed1V1bXUNEiqkbyVlNSuV98/fkqK/9pJTUrxnrGO1x4oNbX6ctucnAzmTWAnTJgAZAIbQx4fSxOYwLQd0/gs7zPKXBk7OuzgxrE3xuz8v1j/C3958jcmM77v+Jid2znHHavvIO9QHoUVhWQPzua4Hsc1fGAE5eZuZsKE67wLpa2Zwe8W/I5lO5fVWz/Fk0LvzH70yehHt7SudEhqT4eUbLJTsslIziDVUkn1pJHmSSHNUkn1pGJt2kG79pgZ/v/yD+ApOYw5SDIDZ96RPQF1qlRmtse1C34P8uzdFdZ7EEBldkdcevB7kGd3HtbA/0f+4zvm4Gp8fiTt3Bb2CsQVnbtCSsB7kKskaWf4/bIV3XpC4L0r5WUk7Q5v6A5mVHTrFbyr5DCefaGH7qxauZKhxx6LS0qiskvw2hhWXIjnQHi/SbvUNCo7dQk+/tBBPAXh9em7NulUdqiRuxzcH5S7hDw+ox2VWR2C9nn278EOh7egnK694Gtv/G/7h31sTfFM9BcBA82sH95fRK6gegx9lVnAtcDHwCXAf51zzsxmAS+Y2cN4b8YdCHyKNw+s1abvmHm+Nqb72nwj1DkaCj6tbRtOOuuUJr/45glvzuz69Wq4SkihF+Zp0NHhTwdYp/5dm3c8zTx/Pcdv27uNAcMHNLPtaMtquEpI9a+8Gp7m3ljY3Bs20xquElJKw1VakGtHXutfHfblVS/HLNF3zsVlas0qZsbEfhN5YfkLAMzbMi/miT54bwS+Zc4tQWsoVOmV1YtJR0/i5KNO5qReJ3F0x6NJ9sTzI7tKv/geP6KZx488unnHM6jhKiGFvh+lLK2Skyec1sxzREuc/+3jfXzcr72midv0mr6e9VuAd4DVwMvOuZVmdreZfctX7Qmgk+9m29uAO3zHrgRexnvj7r+BHzrnKupr09fWL4HbfG118rVd7zlERFq7S4Ze4u/Bzd2Sy9eHvo7JefMO5VFQ6v1lLLtNNl0yujRwRORN7DvRX563ZV5Mz11aUcqTm5/khH+dEJTkpyen8/3R3+eTqZ/w5a1f8s9v/ZPrRl3H4M6DEyTJF5GWJq7vHM65t4G3a+z7dUD5MHBpPcfeC9wbTpu+/ZuonpkncH+95xARac16ZPbg1D6nsmDrAipdJa+tfo2bj7856uet2Ztf7/CpKApM9N/f+j7lleUxSaZ3HtrJRS9dxMfbPvbvS01K5Zbjb+FXp/6Kzm07hzhaRKRxtDKuiMgR7LKhl/nLL698OSbnjOewfja2fwAAEplJREFUnSr9O/Snd5Z3GGRBaQFLdiyJ+jmX5i1l7ONjg5L80/qcxsqbV/LQOQ8pyReRiFOiLyJyBLt46MX+xbIWbF1AXkHjJqFrinhOrVmlapx+lWgP3/n4q4+Z+PREthd4b8jz4OH+M+9n3rXzGNAx0e/tEZGWSom+iMgRrFu7bozv453xxuF4ffXrUT9nIvToA5ze93R/ee6muVE7z4KtCzjr2bPIL/HOeNI+rT1/GP4HfnHyL2K+IrGIHFn0DiMicoS7+JiL/eU3170Z9fMlSqJ/Zv8z/eUFWxdwsORgiNpNszRvKee/cD6FZd7Fi3La5jD/uvmM61jrljERkYhToi8icoQ7f9D5/vK8LfMoKCkIUbt5DpUe4quD3jUikz3J9O/Q9Pmhm6tnVk9GdxsNQFllGe9seCei7W/ct5Fznz/XP8NQt3bdmH/dfEZ2GxnR84iI1EeJvojIEa5Pdh9GdB0BeKd+fHfTu1E717q96/zlAR0HkJIU37UJLhh0gb/81vq3ItbuwZKDnP/i+Xxd6J2yNLtNNv/5zn84Jif0POoiIpGkRF9ERIIS3mgO30mUYTtVLhhc/bpnr5tNeWV4K2eGUukqmTJjiv+1tkluw5tXvsnwrsOb3baISGMo0RcRkaBEf/a62VRUVkTlPIkw406gMd3H0CPTu9r33uK9zNvc/Nl37p5/N7PWzvJvP/GtJzjlqHitpC4iRzIl+iIiwvE9j6drRlcAdhft5tPtn0blPInWo+8xT9BaAs8tf65Z7c1cM5P/nf+//u2fnfQzrhp+VbPaFBFpKiX6IiKCxzxBN+XOXj87KudZtXuVv5wIiT7A1SOu9pdfX/06RWVFTWpn1e5VTJkxxb99Zv8zue/M+5odn4hIUynRFxERAL458Jv+8tvr3454+6UVpazdu9a/PTRnaMTP0RTHdT+OQZ0GAd5ZgWaumdnoNvYX72fy9MkcKj0EQL/sfky/eDrJnuSIxioi0hhK9EVEBPD2QKd4vLPgLN25lB0FOyLa/po9a/w3u/bL7kdmWmZE228qM2PKiOqe+L8s/Eujjq+orODK165kw74NALRNacvMK2bSqW2niMYpItJYSvRFRASAzLRMTutzmn97zvo5EW1/xa4V/nKizUBz/ZjrSU1KBWDh9oV8/NXHYR97x9w7eGdj9Rz8/zf5//zTlYqIxJMSfRER8Ttv4Hn+cqTH6S//erm/PCxnWETbbq6u7bpy9fDqsfr3fRje2PrnvniOBz9+0L9916l3cdmxl4U4QkQkdpToi4iIX+A4/Xc3vUtpRWnE2l6+qzrRT7QefYCfnvhTf3nW2ln8d/N/Q9afu2kuU2dN9W9fMOgC7p54d9TiExFpLCX6IiLiN6jTIPp36A94b0z94MsPItZ2UKLfJfES/eFdh3PNyGv82z+e82OKy4rrrPvBlx8wefpk/xehoTlDee6i5/CYPlZFJHHoHUlERPzMjPMGBAzfWReZ4Tv5h/P5Mv9LAFI8Kf5ZbhLNH874AxkpGQCs3L2SqbOm4pwLqjNzzUwmPTfJPw1n76zezL5qNllpWTGPV0QkFCX6IiIS5JuDAqbZ3BCZaTZX7l7pLw/pPISUpJSItBtpPTJ78NDZD/m3X1zxIpe/ejmf7/ycj776iGtmXMOFL11IYVkhAF0zujL3mrn0ze4bp4hFROqnCX5FRCTI+D7jSU9Op7i8mDV71rBp/yb/cJ6mCrwRNxHH5we6ceyNLNu5jMeWPAbAK6te4ZVVr9Sq179Df9668q2E/XVCREQ9+iIiEiQ9JZ3T+53u347E4lmJPj6/pr+e91d+cNwP6n3+omMu4tPvf8oxOcfEMCoRkcZRoi8iIrVEepXcwER/WJfEmlqzLsmeZKadP433rnmPi4+5mH7Z/RjUaRDfG/U9cq/N5bXLXtOCWCKS8DR0R0REagmcT3/elnkUlRXRNqVtk9qqdJUszVvq3x7ZdWSz44uV0/udHvTrhohIS6IefRERqaVPdh+OzTkWgMPlh5m3eV6T21q/dz0FpQUAdMnoQq+sXhGJUUREQotLom9mHc3sXTNb7/uzQz31rvXVWW9m1wbsP87MlpvZBjP7i5lZqHbN7Goz+8L3+MjMRga0tcXX1jIzWxzt1y4i0lJEapXcJXlL/OXjuh+H7y1bRESiLF49+ncA7znnBgLv+baDmFlH4DfACcA44DcBXwimATcAA32PSQ20uxkY75wbAfwOeLzG6SY650Y558ZG6PWJiLR4gYn+2+vfrjWffLiW7AhO9EVEJDbilehPBp72lZ8Gvl1HnXOAd51z+5xz+4F3gUlm1h3Ics597LyfOs8EHF9nu865j3xtAHwC6HdjEZEGnNz7ZP8iUFvzt7J6z+omtbM4r/rH0uN6KNEXEYmVeCX6XZ1zeQC+P7vUUacn8FXA9jbfvp6+cs394bY7FZgTsO2A/5jZEjO7oQmvRUSkVUpJSuHso8/2bzdlldyaN+KO7aEfTkVEYiVqs+6Y2VygWx1P3RVuE3XscyH2hxPTRLyJ/ikBu092zu0wsy7Au2a2xjm3oJ7jb8A7ZIicnBxyc3PDOa0cIQ4dOqRrQmpp6ddF/4rqhbJeWPQCx5cd36jjtxRu8d+I2yGlA+uXrGeDbYhojC1RS78uJDp0XUikRS3Rd86dWd9zZva1mXV3zuX5huLsqqPaNmBCwHYvINe3v1eN/Tt85XrbNbMRwL+Ac51zewPi3OH7c5eZzcB7P0Cdib5z7nF84/sHDx7sJkyYUFc1OULl5uaia0JqaunXxZBDQ3hg7QMArChYwegTR9O+Tfuwj//nkn/6y6f1P42JEydGPMaWqKVfFxIdui4k0uI1dGcWUDWLzrXAG3XUeQc428w6+G7CPRt4xzckp8DMTvTNtnNNwPF1tmtmRwGvA1Occ+uqTmBmGWaWWVX2nWNF5F6miEjL1q1dN/8NtOWV5by76d1GHf/BVx/4y6ccdUqImiIiEmnxSvTvA84ys/XAWb5tzGysmf0LwDm3D+8MOYt8j7t9+wBuwts7vwHYSPWY+zrbBX4NdAIerTGNZlfgAzP7HPgUmO2c+3eUXrOISIvUnFVyP/zyQ39Zib6ISGzFZWVc39CZM+rYvxj4fsD2k8CT9dSrtYZ6iHa/H9huwP5NQMtZolFEJA7OG3gedy+4G/Am+pWuEo813E+UV5DHxv0bAWiT3IYx3cdENU4REQmmlXFFRCSksT3G0rltZwC+Lvw6aBadUD78qro3f1zPcaQmpUYlPhERqZsSfRERCSnJk8S5A871b4e7Su7cTXP95VN6a9iOiEisKdEXEZEGBa6SO3PNzAbrO+eYs6F6yZJzBpwTlbhERKR+SvRFRKRBkwZMIi0pDYClO5eyaveqkPVX7V7Fl/lfApCVlsVJvU6KeowiIhJMib6IiDQou002Fwy+wL/97OfPhqwf2Jt/Vv+zSElKiVpsIiJSNyX6IiISlmtGXOMvP/vFs1RUVtRb9611b/nLgcN+REQkdpToi4hIWCYNmOSffWd7wfagZD7Q1gNbmb91PgAe8wTdyCsiIrGjRF9ERMKSkpTC1NFT/dsPf/JwnfWe+fwZf/ms/mfRPbN71GMTEZHalOiLiEjYbhl3C8ke71qLC7YuYPGOxUHPO+d46vOn/NvfHfXdWIYnIiIBlOiLiEjYemX14tKhl/q3b//P7Tjn/Nuz1s5i0/5NALRPa8/kIZNjHqOIiHgp0RcRkUa569S7SLIkAOZvnc9LK18CoLSilDv/e6e/3tTRU2mT3CYuMYqIiBJ9ERFppGO7HMuPT/ixf/sHb/2Af2/4Nz+c/UP//PoZKRn8/OSfxytEEREBkuMdgIiItDy/Gf8bXl31Kl8d/Ir8knzOfT54Zp17Tr+Hbu26xSk6EREB9eiLiEgTtG/TnllXzqJjesdaz101/KqgHn8REYkPJfoiItIko7qNYtH1i7hi2BV0zejKkM5DeOjsh3jm28/gMX28iIjEm4buiIhIk/Xv0J8XL34x3mGIiEgd1OUiIiIiItIKKdEXEREREWmFlOiLiIiIiLRCSvRFRERERFohJfoiIiIiIq2QEn0RERERkVZIib6IiIiISCukRF9EREREpBVSoi8iIiIi0gop0RcRERERaYXMORfvGFokMysA1sY7DkkonYE98Q5CEo6uC6mLrgupi64Lqctg51xmUw5MjnQkR5C1zrmx8Q5CEoeZLdY1ITXpupC66LqQuui6kLqY2eKmHquhOyIiIiIirZASfRERERGRVkiJftM9Hu8AJOHompC66LqQuui6kLroupC6NPm60M24IiIiIiKtkHr0RURERERaISX6IZjZJDNba2YbzOyOOp5PM7OXfM8vNLO+sY9SYi2M6+I2M1tlZl+Y2Xtm1icecUpsNXRdBNS7xMycmWlmjSNAONeFmV3me89YaWYvxDpGib0wPkeOMrN5ZrbU91lyXjzilNgxsyfNbJeZrajneTOzv/iumS/MbEw47SrRr4eZJQF/B84FhgJXmtnQGtWmAvudcwOAPwH3xzZKibUwr4ulwFjn3AjgVeCB2EYpsRbmdYGZZQI/BhbGNkKJh3CuCzMbCPwKONk5dyxwa8wDlZgK8/3i/wEvO+dGA1cAj8Y2SomDp4BJIZ4/Fxjoe9wATAunUSX69RsHbHDObXLOlQLTgck16kwGnvaVXwXOMDOLYYwSew1eF865ec65It/mJ0CvGMcosRfO+wXA7/B+8Tscy+AkbsK5Lq4H/u6c2w/gnNsV4xgl9sK5LhyQ5Su3B3bEMD6JA+fcAmBfiCqTgWec1ydAtpl1b6hdJfr16wl8FbC9zbevzjrOuXIgH+gUk+gkXsK5LgJNBeZENSJJBA1eF2Y2GujtnHsrloFJXIXzfjEIGGRmH5rZJ2YWqkdPWodwrovfAt8xs23A28CPYhOaJLDG5h+AVsYNpa6e+ZpTFIVTR1qXsP/Nzew7wFhgfFQjkkQQ8rowMw/e4X3XxSogSQjhvF8k4/0pfgLeX//eN7NhzrkDUY5N4iec6+JK4Cnn3ENmdhLwrO+6qIx+eJKgmpRzqke/ftuA3gHbvaj905m/jpkl4/15LdTPLtLyhXNdYGZnAncB33LOlcQoNomfhq6LTGAYkGtmW4ATgVm6IbfVC/dz5A3nXJlzbjOwFm/iL61XONfFVOBlAOfcx0AboHNMopNEFVb+UZMS/fotAgaaWT8zS8V7M8ysGnVmAdf6ypcA/3VamKC1a/C68A3R+AfeJF/jbY8MIa8L51y+c66zc66vc64v3ns3vuWcWxyfcCVGwvkcmQlMBDCzzniH8myKaZQSa+FcF18CZwCY2TF4E/3dMY1SEs0s4Brf7DsnAvnOubyGDtLQnXo458rN7BbgHSAJeNI5t9LM7gYWO+dmAU/g/TltA96e/CviF7HEQpjXxR+BdsArvnuzv3TOfStuQUvUhXldyBEmzOviHeBsM1sFVAA/d87tjV/UEm1hXhc/A/5pZj/FOzzjOnUktm5m9iLeIXydffdm/AZIAXDOPYb3Xo3zgA1AEfDdsNrVdSMiIiIi0vpo6I6IiIiISCukRF9EREREpBVSoi8iIiIi0gop0RcRERERaYWU6IuIiIiItEJK9EVEREREWiEl+iIiUicz62Rmy3yPnWa2PWD7oyidc7SZ/SvE8zlm9u9onFtEpLXRglkiIlIn38JNowDM7LfAIefcg1E+7Z3APSFi2m1meWZ2snPuwyjHIiLSoqlHX0REGs3MDvn+nGBm883sZTNbZ2b3mdnVZvapmS03s6N99XLM7DUzW+R7nFxHm5nACOfc577t8QG/ICz1PQ8wE7g6Ri9VRKTFUqIvIiLNNRL4CTAcmAIMcs6NA/4F/MhX5xHgT86544GLfc/VNBZYEbB9O/BD59wo4FSg2Ld/sW9bRERC0NAdERFprkXOuTwAM9sI/Me3fzkw0Vc+ExhqZlXHZJlZpnOuIKCd7sDugO0PgYfN7HngdefcNt/+XUCPyL8MEZHWRYm+iIg0V0lAuTJgu5LqzxkPcJJzrpj6FQNtqjacc/eZ2WzgPOATMzvTObfGVydUOyIigobuiIhIbPwHuKVqw8xG1VFnNTAgoM7Rzrnlzrn78Q7XGeJ7ahDBQ3xERKQOSvRFRCQWfgyMNbMvzGwV8IOaFXy99e0Dbrq91cxWmNnneHvw5/j2TwRmxyJoEZGWzJxz8Y5BREQEADP7KVDgnAs1l/4CYLJzbn/sIhMRaXnUoy8iIolkGsFj/oOYWQ7wsJJ8EZGGqUdfRERERKQVUo++iIiIiEgrpERfRERERKQVUqIvIiIiItIKKdEXEREREWmFlOiLiIiIiLRC/x8EruHIUK6CWAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E_3pt =  9.693346990256778e-08\n"
     ]
    }
   ],
   "source": [
    "# Compare FD Seismogram with analytical solution\n",
    "# ---------------------------------------------- \n",
    "# Define figure size\n",
    "rcParams['figure.figsize'] = 12, 5\n",
    "\n",
    "# plot solution of 3-point FD operator\n",
    "plt.plot(time, seis_3pt, 'b-',lw=3,label=\"FD solution (3pt operator)\")\n",
    "\n",
    "# PLOT SOLUTION OF 5-POINT FD OPERATOR HERE!\n",
    "\n",
    "# plot analytical solution\n",
    "Analy_seis = plt.plot(time,seis_ana,'r--',lw=3,label=\"Analytical solution\") \n",
    "\n",
    "# plot difference between analytical and 3pt operator FD solution\n",
    "plt.plot(time, 10*(seis_3pt-seis_ana), 'g-',lw=3,label=\"3pt operator error (10x)\")\n",
    "\n",
    "# PLOT DIFFERENCE BETWEEN ANALYTICAL AND 5PT OPERATOR SOLUTION HERE!\n",
    "\n",
    "plt.xlim(time[0], time[-1])\n",
    "plt.title('Seismogram comparison')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# Calculate L2-norm of the data residuals between FD and analytical solution\n",
    "# 3-point operator\n",
    "E_3pt = np.sum((seis_3pt-seis_ana)**2)\n",
    "print(\"E_3pt = \", E_3pt)\n",
    "\n",
    "# ADD L2-NORM for 5-POINT OPERATOR HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "157 ms ± 614 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "dx = 2.0 # grid point distance in x-direction\n",
    "dt = 0.0010  # time step\n",
    "\n",
    "# Benchmark runtime of 3pt-FD code\n",
    "op = 3 # use 3-point operator\n",
    "%timeit time, seis_3pt, seis_ana = FD_1D_acoustic(dt,dx,op)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BENCHMARK RUNTIME OF 5PT-FD CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAusAAAFNCAYAAAC5VCFTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmcFNW5//HPw7AvgiyCLLIzyCIgAtG4oETFGyOuV4waJSYkMSS58RcTTPLzGiKJJrlqcqMxJrjcxADq/akYMUajY+LGIoKyiA6LMiKyK9uwDM/vj6qhe3q6Z3qge6p7+vt+vfo1p5Zz6qmemplnTp86Ze6OiIiIiIjknkZRByAiIiIiIskpWRcRERERyVFK1kVEREREcpSSdRERERGRHKVkXUREREQkRylZFxERERHJUUrWRUQywMx2mlmfqOOQujGzZWY2Nuo4RERSMc2zLiISMLNTgV8Ag4EKYAXwH+6+INLARESkYDWOOgARkVxgZkcBfwW+ATwCNAVOA/ZGGVemmFmRu1dEHUeuMLPG7n4g6jhERGqjYTAiIoEBAO4+090r3H2Pu//d3d+q3MHMvmxmK8xsm5k9a2Y947a5mfULy/9mZsvNbIeZfWhm3wvXjzWzMjP7vpltNLOPzOzCcP93zWyrmf0wrs1mZnaXma0PX3eZWbO47d8P21hvZl9JiOFBM/udmc01s13AmWb2eTN708w+NbN1ZnZLXFu9wvqTwm3bzOzrZjbKzN4ys+1m9ttUb56ZFZnZD81sVXjeb5hZj3DbKWa2wMw+Cb+eElevxMxuNbNXw6FET5lZBzN7OIxzgZn1Snifv21mq81ss5n90swahdv6mtkLZrYl3PawmbWLq7vWzH5gZm8Bu8yscbjuc+H20Wa2MDzux2Z2R1zdC8IhM9vDmI9PaPd74fv0iZnNNrPmNV5tIiJpUrIuIhJ4F6gws4fM7DwzOzp+o5ldCPwQuBjoBPwLmJmirRnA19y9DTAEeCFuWxegOdANuBn4A3AVMJKgJ//muLHvPwI+AwwHhgGjgR+H8YwHbgA+B/QDzkgSxxeB6UAb4GVgF/AloB3weeAb4XnFGwP0By4H7gpj+BzB0KB/N7NkxyGM5Qrg34CjgC8Du82sPfA08BugA3AH8LSZdYirOxG4OnxP+gKvAQ8A7QmGIv1nwrEuAk4CTgQmhMcCMODnQFfgeKAHcEtC3SvCc2+XpGf918Cv3f2oMI5HAMxsAMH3+j8IvvdzgafMrGlc3X8HxgO9gROAa1O8TyIidaJkXUQEcPdPgVMBJ0igN5nZHDPrHO7yNeDn7r4iTPJ+BgyP712Psx8YZGZHufs2d1+UsG26u+8HZgEdCRLEHe6+DFhGkOwBXAlMc/eN7r4J+AlBUgtBcviAuy9z993htkRPuvsr7n7Q3cvdvcTd3w6X3yJIQBOT75+G+/6dILmfGR7/Q4J/UEakeAu/AvzY3Vd6YIm7byFIjN9z9z+5+wF3nwm8A3whru4D7r7K3T8BngFWufvz4fv8aJJj3u7uW939A4J/KK4AcPdSd3/O3feG79cdSc7vN+6+zt33JDmH/UA/M+vo7jvd/fVw/eXA02Hb+4FfAS2AU+Lq/sbd17v7VuApgn+wRESOmJJ1EZFQmIhf6+7dCXrEuxIkgwA9gV+HwyC2A1sJenK7JWnqEoIe5vfN7CUzOzlu25a4seOVCePHcdv3AK3Dclfg/bht74frKreti9sWX066zszGmNmLZrbJzD4Bvk7wz0K8xFhSxZaoB7AqyfrEcyBcjn/f6nrM+PM69J6Y2TFmNsuCoUefAn+m+vkle58qXUcwHOqdcPjN+cnOwd0Phu3En8OGuPLuJDGLiBwWJesiIkm4+zvAgwRJOwTJ2dfcvV3cq4W7v5qk7gJ3nwAcAzxBOJziMKwn+Ceh0nHhOoCPgO5x23okO42E5b8Ac4Ae7t4WuJfgH45MWEcwdCRR4jlAcB4fHsGx4s81/j35OcE5nxAOZbmK6ueXcgo0d3/P3a8g+L7dDjxmZq1IOAczszCGIzkHEZG0KFkXEQHMbKCZ/R8z6x4u9yAYXlE5FOJe4CYzGxxub2tmlyVpp6mZXWlmbcMhE58STAN5OGYCPzazTmbWkWCM+5/DbY8Ak8zseDNrGW6rTRtgq7uXm9logjHtmfJH4Kdm1t8CJ4Tj0ucCA8zsi+ENnZcDgwhm3jlcN5rZ0eH36DvA7HB9G2AnsN3MugE31qVRM7vKzDqFPefbw9UVBO/1581snJk1Af4PwSxB1f5RExHJNCXrIiKBHQQ3V86zYPaU14GlBIkZ7v44QW/rrHCIxVLgvBRtXQ2sDff7OkEP7+G4FVgIvAW8DSwK1+HuzxDctPkiUEpwUybUPNXk9cA0M9tBkNwfbo9/MneE7f2d4B+UGUCLcNz6+QTv4xbg+8D57r75CI71JPAGsJjg5tUZ4fqfENx0+km4/v/Vsd3xwDIz20lws+nEcPz+SoLv4X8DmwnG23/B3fcdwTmIiKRFD0USEWkAwqkElwLNGvL84WbmQH93L406FhGR+qCedRGRPGVmF4XDbo4m6PV/qiEn6iIihUjJuohI/voasIlgFpYKgqeviohIA6JhMCIiIiIiOUo96yIiIiIiOUrJuoiIiIhIjmocdQBRateunffr1y/qMCTH7Nq1i1atWkUdhuQYXReSjK4LSUbXhSTzxhtvbHb3TnWtV9DJeufOnVm4cGHUYUiOKSkpYezYsVGHITlG14Uko+tCktF1IcmY2fuHU0/DYEREREREcpSSdRERERGRHKVkXUREREQkRxX0mHURERGRTDMz1qxZQ3l5edShSASaN29O9+7dadKkSUbaU7IuIiIikkGtWrWiTZs29OrVCzOLOhypR+7Oli1bKCsro3fv3hlpU8NgRERERDKoqKiIDh06KFEvQGZGhw4dMvqpipJ1ERERkQxTol64Mv29V7IuIiIi0gA9/vjjmBnvvPPOEbVz7bXX8thjj9W4z89+9rMqy6eccsphHeuWW27hV7/61WHVrVRSUsL5559f4z7bt2/nnnvuObS8fv16Lr300iM6brYoWRcRERFpgGbOnMmpp57KrFmzsn6sxGT91Vdfzfoxj0Rist61a9da/yGJipJ1ERERkQZm586dvPLKK8yYMaNKsl75dNVLL72UgQMHcuWVV+LuAEybNo1Ro0YxZMgQJk+efGh9pX/84x9cdNFFh5afe+45Lr74YqZOncqePXsYPnw4V155JQCtW7c+tN8vfvELhg4dyrBhw5g6dSoAf/jDHxg1ahTDhg3jkksuYffu3TWez6OPPsqQIUMYNmwYp59+OgDl5eVMmjSJoUOHMmLECF588cVq9RJ76ocMGcLatWuZOnUqq1atYvjw4dx4442sXbuWIUOG1Njugw8+yMUXX8z48ePp378/3//+92v5LmSGknUREcmMVatg7FjYsiXqSEQK3hNPPMH48eMZMGAA7du3Z9GiRYe2vfnmm9x1110sX76c1atX88orrwAwZcoUFixYwNKlS9mzZw9//etfq7R51llnsWLFCjZt2gTAAw88wKRJk7jtttto0aIFixcv5uGHH65S55lnnuGJJ55g3rx5LFmy5FCCe/HFF7NgwQKWLFnC8ccfz4wZM2o8n2nTpvHss8+yZMkS5syZA8Ddd98NwNtvv83MmTO55ppr0r6x87bbbqNv374sXryYX/7yl1W21dTu4sWLmT17Nm+//TazZ89m3bp1aR3vSChZFxGRunvhBaioiC2vXg1nngkvvQTf+lZ0cYnkoltuAbP0XpMnV68/eXLVfW65pdZDzpw5k4kTJwIwceJEZs6ceWjb6NGj6d69O40aNWL48OGsXbsWgBdffJExY8YwdOhQXnjhBZYtW1alTTPj6quv5s9//jPbt2/ntdde47zzzqsxjueff55JkybRsmVLANq3bw/A0qVLOe200xg6dCgPP/xwtWMl+uxnP8u1117LH/7wByrC3z0vv/wyV199NQADBw6kZ8+evPvuu7W+N7Wpqd1x48bRtm1bmjdvzqBBg3j//feP+Hi10TzrIiJSNwsXwrhxMGwY3HVX0Jv+9ttQ2cM0cyZ897swalSkYYoUqi1btvDCCy+wdOlSzIyKigrMjF/84hcANGvW7NC+RUVFHDhwgPLycq6//noWLlxIjx49uOWWW5L2Uk+aNIkvfOELNG/enMsuu4zGjWtOJd096ewo1157LU888QTDhg3jwQcfpKSkpMZ27r33XubNm8fTTz/N8OHDWbx4cbVhOsk0btyYgwcPHlpOp+e9pnaTvXfZpp51ERGpm7vuCr4uWQK//W1QnjABLrssts/06fUfl4gA8Nhjj/GlL32J999/n7Vr17Ju3Tp69+7Nyy+/nLJOZRLbsWNHdu7cmfJmy65du9K1a1duvfVWrr322kPrmzRpwv79+6vtf84553D//fcfGpO+detWAHbs2MGxxx7L/v37qw2dSWbVqlWMGTOGadOm0bFjR9atW8fpp59+qO67777LBx98QHFxcZV6vXr1OjQEaNGiRaxZswaANm3asGPHjqTHSqfd+qRkXURE0rdpEzz6aGz5hz+MleM/mn/qKaiHj4dF8sItt4B7eq/77qte/777qu5TyzCYmTNnVrkRFOCSSy7hL3/5S8o67dq146tf/SpDhw7lwgsvZFQNn4xdeeWV9OjRg0GDBh1aN3nyZE444YRDN5hWGj9+PBdccAEnnXQSw4cPP3Sz509/+lPGjBnD2WefzcCBA2s8H4Abb7yRoUOHMmTIEE4//XSGDRvG9ddfT0VFBUOHDuXyyy/nwQcfrNLzXXneW7duZfjw4fzud79jwIABAHTo0IHPfvazDBkyhBtvvLFKnXTarU+WzkcIDVVxcbGvXLky6jAkx1TeKS8ST9dF6PbbIZzNgTFj4PXXq24/5xx47rmgfNtt8IMf1G989UzXhSTz5ptvMmLEiKjDyJopU6YwYsQIrrvuuqhDyVkrVqzg+OOPr7LOzN5w95Pq2pZ61kVEJH3xH1d/4xvVt191Vaz8yCPZj0dE6tXIkSN56623uCr+Z12ySjeYiohIetasCW4kBWjWDC65pPo+F1wATZvCvn2waFEwnWPfvvUbp4hkzRtvvBF1CAVHPesiIpKep56KlceNg7iHnhzSrl0wFKZS/Ph2ERGpMyXrIiKSniefjJUvuCD1fvGzwjzzTPbiEclhhXxPYKHL9PdeybqIiNTu00/hn/+MLZ9/fup943vWFy8OhsSIFJCKigq2bNmihL0AuTtbtmyhefPmGWtTY9ZFRKR2L78MlQ//GD4cunVLvW+XLnDHHcFDk045JRjDLlJAdu3axY4dO9i0aVPUoUgEmjdvTvfu3TPWnpJ1ERGpXdeuwewvL70EZ55Z+/7f/W72YxLJUe5O7969ow5DGoisDoMxs/FmttLMSs1sapLtzcxsdrh9npn1itt2U7h+pZmdG7f+fjPbaGZLUxzze2bmZtYxG+ckIlKQhg+He+6BZcsgfKiJiIhkX9aSdTMrAu4GzgMGAVeY2aCE3a4Dtrl7P+BO4Paw7iBgIjAYGA/cE7YH8GC4LtkxewBnAx9k9GRERCSmkW53EhGpL9n8jTsaKHX31e6+D5gFTEjYZwLwUFh+DBhnZhaun+Xue919DVAatoe7/xPYmuKYdwLfB3RHh4hILti2DZ59NphvXURE6iybyXo3YF3cclm4Luk+7n4A+ATokGbdKszsAuBDd19yZGGLiEhG3HwztG8P48fD7NlRRyMikpeyeYOpJVmX2OOdap906sYaMWsJ/Ag4J9U+cftOBiYDdOrUiZKSktqqSIHZuXOnrgupppCvi4G33UaT7dvZMXAgH33+8+zt1Cmtel327WNgWN789NMsPeWU7AUZkUK+LiQ1XReSSdlM1suAHnHL3YH1KfYpM7PGQFuCIS7p1I3XF+gNLAlG0dAdWGRmo919Q/yO7n4fcB9AcXGxjx07tm5nJQ1eSUkJui4kUcFeF+5w+eWwcSMd5s2j1w9+AMcfn17dLl3g9tsB6Pjee4w94wywZH0x+atgrwupka4LyaRsDoNZAPQ3s95m1pTghtE5CfvMAa4Jy5cCL3jwBIE5wMRwtpjeQH9gfqoDufvb7n6Mu/dy914Eyf6JiYm6iIjU0bp1sHFjUG7TBoqL0687YAC0axeUN22CNWsyH5+ISAOXtWQ9HIM+BXgWWAE84u7LzGxaOL4cYAbQwcxKgRuAqWHdZcAjwHLgb8A33b0CwMxmAq8BxWZWZmbXZescREQK3oIFsfLIkXWbCaZRIxgzJrY8P2Wfi4iIpJDVhyK5+1xgbsK6m+PK5cBlKepOB6YnWX9FGsftVddYRUQkifhkfdSoutcfOTKYDQZg8WKYODEzcYmIFAhNlisiIqktXBgrH06yPnx4rPzmm0cej4hIgVGyLiIiybnDkrjZcE88se5tjBgRK7/5ZtCmiIikTcm6iIgkt2EDbN4clFu3ht69695Gnz5BXQhuMv3oo8zFJyJSAJSsi4hIcvG96iecULebSys1agTDhsWWFy8+8rhERApIVm8wFRGRPPbWW7HyCSccfjsjRwa96sOHx6ZyFBGRtChZFxGR5OJ71uN7x+vqrrsa3MOQRETqi5J1ERFJ7ve/hylTgh72ceMOvx0l6iIih03JuoiIJNe6NZx8cvASEZFI6AZTEREREZEcpZ51ERHJvnnzoKQEli+Ha66Bs86KOiIRkbygZF1ERKpbtw5atICOHTPT3iOPwB13BOW+fZWsi4ikScNgRESkuhtvhE6dgtdf/3rk7Q0aFCsvX37k7YmIFAgl6yIiUt2KFcHXzZuhQ4cjb2/w4FhZybqISNqUrIuISFUVFbByZWz5+OOPvM34NlauhAMHjrxNEZECoGRdRESqWrsW9u4Nyl26ZOapo23bQrduQXnfPli16sjbFBEpAErWRUSkqvhhKpnoVa+koTAiInWmZF1ERKqqHK8OmU3WdZOpiEidKVkXEZGqlKyLiOQMJesiIlJVtpL1gQNj5XffzVy7IiINmJJ1ERGJcc9esj5gQKz87rvBsUREpEZ6gqmIiMRs2ACffhqUjzoKjj02c20fcwxcfz307h0k7u5glrn2RUQaICXrIiISs3kz9O0bTN84YEBmk2kzuPvuzLUnIlIAlKyLiEjM0KFQWhrMhb5tW9TRiIgUPI1ZFxGR6po2hc6do45CRKTgZTVZN7PxZrbSzErNbGqS7c3MbHa4fZ6Z9YrbdlO4fqWZnRu3/n4z22hmSxPa+qWZvWNmb5nZ42aWgUfuiYhI1ugGUxGRWmUtWTezIuBu4DxgEHCFmQ1K2O06YJu79wPuBG4P6w4CJgKDgfHAPWF7AA+G6xI9Bwxx9xOAd4GbMnpCIiJy5NatgwsuCKZxHDUq6mhERHJeNnvWRwOl7r7a3fcBs4AJCftMAB4Ky48B48zMwvWz3H2vu68BSsP2cPd/AlsTD+buf3f3A+Hi60D3TJ+QiEiD5g6/+Q08/TSsXJmdnu9WreCpp4L2ly2DgwczfwwRkQYkmzeYdgPWxS2XAWNS7ePuB8zsE6BDuP71hLrd6nDsLwOz6xqwiEhBW78evvOdoHz00bC1Wr/IkWvfHjp2DGadKS+HsjI47rjMH0dEpIHIZrKebL6vxG6aVPukUzf5Qc1+BBwAHk6xfTIwGaBTp06UlJSk06wUkJ07d+q6kGoK4bpot3gxw8Pyp507syhL5zuic2fabt4MwOJHH2X7yJFZOU59KITrQupO14VkUjaT9TKgR9xyd2B9in3KzKwx0JZgiEs6dasxs2uA84Fx7sk/v3X3+4D7AIqLi33s2LHpnIsUkJKSEnRdSKKCuC5KSw8Vjxo5Mnvne9JJwRAYYHjLlpDH72tBXBdSZ7ouJJOyOWZ9AdDfzHqbWVOCG0bnJOwzB7gmLF8KvBAm2XOAieFsMb2B/sD8mg5mZuOBHwAXuPvuDJ6HiEhheO+9WLl//+wdJ77t+GOKiEg1WUvWw5s9pwDPAiuAR9x9mZlNM7MLwt1mAB3MrBS4AZga1l0GPAIsB/4GfNPdKwDMbCbwGlBsZmVmdl3Y1m+BNsBzZrbYzO7N1rmJiDRI8Ylzv37ZO46SdRGRtGX1CabuPheYm7Du5rhyOXBZirrTgelJ1l+RYv8s/mURESkA9dWzHv+PwKpV2TuOiEgDoCeYiohIMIVi3Jj1rCbrffrEymvWaPpGEZEaKFkXEZFg2sby8qDcvn0wdWO2tGsXa7+8HD76KHvHEhHJc0rWRUQEVq+Olfv2zf7x4o+hoTAiIilldcy6iIjkiTVrYuX4YSrZ8l//BY0bB8fq3Dn7xxMRyVNK1kVEBHr0gC9+MUjaR4zI/vFOPz37xxARaQCUrIuICJx1VvASEZGcojHrIiIiIiI5Ssm6iIhEZ9s2ePPNqKMQEclZGgYjIiL1zx26dIGNG4PlHTugdetoYxIRyUHqWRcRKXQffBDcXPqjH8Ejj9TPMc2gbdvYcvxsNCIicoiSdRGRQrdiBcycCT/7Gdx9d/0dN36KSM21LiKSlJJ1EZFCF9+r3bt3/R03/sFI8Q9lEhGRQ5Ssi4gUuqiSdfWsi4jUSsm6iEihU8+6iEjOUrIuIlLo1LMuIpKzlKyLiBS6qJL1+GOtXQsVFfV3bBGRPKFkXUSkkO3YAVu2BOWmTaFr1/o7dps2cMwxQXn/figrq79ji4jkCSXrIiKFLL5XvWdPaFTPfxbih8Jo3LqISDVK1kVECllUQ2Aq9ekDRUXBsffsqf/ji4jkuMZRByAiIhGKOlm/5x546CForD9HIiLJ6LejiEghO+ec4Kmlq1fDqafW//Hbtq3/Y4qI5BEl6yIihWzQoOAlIiI5SWPWRURERERylHrWRUQkWu+8EzwUafVq+MpXoEWLqCMSEckZWe1ZN7PxZrbSzErNbGqS7c3MbHa4fZ6Z9YrbdlO4fqWZnRu3/n4z22hmSxPaam9mz5nZe+HXo7N5biIikiHnnx+8vv1tTd8oIpIga8m6mRUBdwPnAYOAK8wscWDkdcA2d+8H3AncHtYdBEwEBgPjgXvC9gAeDNclmgr8w937A/8Il0VEJJWFC4O51ceOhVtuiS6O+LnW42enERGRrPasjwZK3X21u+8DZgETEvaZADwUlh8DxpmZhetnufted18DlIbt4e7/BLYmOV58Ww8BF2byZEREGpxVq+CDD+Cll2Dx4ujiiJ8yUj3rIiJVZDNZ7wasi1suC9cl3cfdDwCfAB3SrJuos7t/FLb1EXDMYUcuIlIIop5jPdmx1bMuIlJFNm8wtSTrPM190ql7WMxsMjAZoFOnTpSUlGSiWWlAdu7cqetCqmmI18WAV16ha1h+78ABPozo/Drt2cPgsLx5/nyW5tH73BCvCzlyui4kk7KZrJcBPeKWuwPrU+xTZmaNgbYEQ1zSqZvoYzM71t0/MrNjgY3JdnL3+4D7AIqLi33s2LHpnY0UjJKSEnRdSKIGeV387GeHiv3POYf+UZ1fq1YwbRoAHXfsyKv3uUFeF3LEdF1IJmVzGMwCoL+Z9TazpgQ3jM5J2GcOcE1YvhR4wd09XD8xnC2mN9AfmF/L8eLbugZ4MgPnICLScOXiMJjVq8Ez8kGqiEiDkLVkPRyDPgV4FlgBPOLuy8xsmpldEO42A+hgZqXADYQzuLj7MuARYDnwN+Cb7l4BYGYzgdeAYjMrM7PrwrZuA842s/eAs8NlERFJpqIC3n8/ttyrV2Sh0KEDtGkTlHftgs2bo4tFRCTHZPWhSO4+F5ibsO7muHI5cFmKutOB6UnWX5Fi/y3AuCOJV0SkYHz4IezfH5Q7dYLWraOLxSzoXX/rrWB5zZogJhERye5DkUREJEflyhCYZDFo+kYRkUOy2rMuIiI5KteS9RNPhI0bgwckde1a+/4iIgVCybqISCHKtWT95puDl4iIVKFkXUSkEH3ve3DRRcGQk379oo5GRERSULIuIlKI2rSB4cODl4iI5CzdYCoiIiIikqPUsy4iIrlh9uxg+sY1a+BXv9KNpiIiKFkXESk8e/bAtm3QpQs0yqEPWH/9a3jttaD81a8qWRcRQcNgREQKzz//Cd26QcuWcNVVUUcTEz8rTfxsNSIiBUzJuohIoalMhPfuhcY59AFrnz6xsh6MJCICKFkXESk8uTbHeiX1rIuIVKNkXUSk0ChZFxHJG0rWRUQKTa4m6xoGIyJSTa3Jupm1NLP/a2Z/CJf7m9n52Q9NRESyIleT9e7dY2PoP/4Ydu+ONh4RkRyQTs/6A8Be4ORwuQy4NWsRiYhI9uzYAVu2BOWmTXNresSiIjjuuNiyhsKIiKSVrPd1918A+wHcfQ9gWY1KRESyIz4B7tkzt+ZZh6pDYZSsi4iklazvM7MWgAOYWV+CnnYREck38QlwfGKcK+KH5WjcuohIWk8w/U/gb0APM3sY+CxwbTaDEhGRLMnV8eqVzjsPjjoq+Edi7NiooxERiVytybq7P2dmi4DPEAx/+Y67b856ZCIiknm7dkGLFrBnT24m6xddFLxERASoIVk3sxMTVn0Ufj3OzI5z90XZC0tERLLiRz+CH/4wmG2ladOooxERkVrU1LP+X+HX5sBJwBKCnvUTgHnAqdkNTUREssIMunSJOgoREUlDyhtM3f1Mdz8TeB840d1PcveRwAigtL4CFBEREREpVOnMBjPQ3d+uXHD3pcDw7IUkIiIF7Sc/gTPPhF69YOHCqKMREYlUOrPBrDCzPwJ/Jpi+8SpgRVajEhGRzFu/Hl59NbixtE8fOProqCNK7s03oaQkKK9aBSedFGk4IiJRSqdnfRKwDPgO8B/A8nBdrcxsvJmtNLNSM5uaZHszM5sdbp9nZr3itt0Url9pZufW1qaZjTOzRWa22MxeNrN+6cQoIlIw/vUvuOyyIPmdlNav8WhornURkUPSmbqxHLgzfKXNzIqAu4GzgTJggZnNcfflcbtdB2xz935mNhG4HbjczAYBE4HBQFfgeTMbENZJ1ebvgAnuvsLMrgd+jOaDFxGJyfU51ivpKaYiIofUmqyb2RrCp5fGc/faHn03GihrJosMAAAgAElEQVR199VhO7OACQQ985UmALeE5ceA35qZhetnufteYI2ZlYbtUUObDhwV7tMWWF/buYmIFJR8SdbjY1OyLiIFLp0x6/GDBZsDlwHt06jXDVgXt1wGjEm1j7sfMLNPgA7h+tcT6nYLy6na/Aow18z2AJ8SPMRJREQqxQ8p6VNbf0uENAxGROSQdIbBbElYdZeZvQzcXEtVS9ZcmvukWp9sjH1lm98F/s3d55nZjcAdBAl81QOaTQYmA3Tq1ImSypuYREI7d+7UdSHVNITrYszy5bQIy/M3b2Z3jp5Po/JyTg/LB99/n3/+4x9QVBRpTKk0hOtCMk/XhWRSOsNg4p9k2oigp71NGm2XAT3ilrtTfWhK5T5lZtaYYPjK1lrqVltvZp2AYe4+L1w/G/hbsqDc/T7gPoDi4mIfO3ZsGqcihaSkpARdF5Io76+LAwdg48ZDi6P//d+hZcsIA6pF587w8cc0qqhgbL9+0LNn1BEllffXhWSFrgvJpHSGwfxXXPkAsAb49zTqLQD6m1lv4EOCG0a/mLDPHOAa4DXgUuAFd3czmwP8xczuILjBtD8wn6DHPVmb24C2ZjbA3d8luAFV00uKiFQqKwsSdgieXprLiToEQ2E+/jgor16ds8m6iEi2pZOsX1d5Q2elMFmuUTgGfQrwLFAE3O/uy8xsGrDQ3ecAM4A/hTeQbiVIvgn3e4TgxtEDwDfdvSI8drU2w/VfBf7XzA4SJO9fTuPcREQKQ76MV6/Upw+8Ht66tGZN8JAkEZEClE6y/hhwYpJ1I2ur6O5zgbkJ626OK5cT3LCarO50YHo6bYbrHwcery0mEZGClC8zwVTSjDAiIkANybqZDSSY57ytmV0ct+kogllhREQkX+Rbz/rll8OJJwZJe9++UUcjIhKZmnrWi4HzgXbAF+LW7wC+ms2gREQkw3r1grPOCnqp+/ePOpraDR0avEREClzKZN3dnwSeNLOT3f21eoxJREQy7atfDV4iIpJXahoG8313/wXwRTO7InG7u387q5GJiIiIiBS4mobBVE59uLA+AhEREanm4EHYsAHatoVWraKORkSk3tU0DOap8OtD9ReOiIhIaNIkmDkT9u6Fp56C88+POiIRkXpX0zCYpwBPtd3dL8hKRCIiklnz5sHs2cHMKqNHw5gxUUeUniZNgkQdNH2jiBSsmobB/KreohARkex59VW4886g/PWv50+yHj/X+urVqfcTEWnAahoG81Jl2cyaAgMJetpXuvu+eohNREQyIb5XOh/mWK8UH6t61kWkQNX6BFMz+zxwL7AKMKC3mX3N3Z/JdnAiIpIB8b3S+fD00kp6iqmISO3JOvBfwJnuXgpgZn2BpwEl6yIi+aAh9KyvXg3uYBZdPCIiEWiUxj4bKxP10GpgY5biERGRTHKv2rOeT8l6hw7QunVQ3rkTtmyJNh4RkQikk6wvM7O5ZnatmV0DPAUsMLOLzeziLMcnIiJHYsMGKC8Pyu3aBa98YaahMCJS8NJJ1psDHwNnAGOBTUB74AuAJr0VEcll+ToEplLiUBgRkQJT65h1d59UH4GIiEgW5OvNpZXUsy4iBS6d2WB6A98CesXvr4ciiYjkgXwdr16pMuYWLWD37mhjERGJQDqzwTwBzCAYq34wu+GIiEhG5fswmKuvhssug86dNROMiBSkdJL1cnf/TdYjERGRzLv8cujRI+hhHzEi6mjqLp9uiBURyYJ0kvVfm9l/An8H9laudPdFWYtKREQyY/z44CUiInkpnWR9KHA1cBaxYTAeLouIiIiISJakk6xfBPRx933ZDkZERKSa7dth1apg/P2oUdCzZ9QRiYjUm3TmWV8CaNCgiIhE4/rr4aSTghtN//GPqKMREalX6fSsdwbeMbMFxMasu7tPyF5YIiJyxO65Bx58EPr1gyuvhM9/PuqIDk/8LDaaa11ECkw6yfp/xpUNOBW4IjvhiIhIxixZAgsWBK/Ro6OO5vDpwUgiUsBqHQbj7i8BnwCfBx4ExgH3ptO4mY03s5VmVmpmU5Nsb2Zms8Pt88ysV9y2m8L1K83s3NratMB0M3vXzFaY2bfTiVFEpMEqLY2V+/ePLo4jFd+zHv+QJxGRApCyZ93MBgATCXrRtwCzAXP3M9Np2MyKgLuBs4EyYIGZzXH35XG7XQdsc/d+ZjYRuB243MwGhcceDHQFng/joYY2rwV6AAPd/aCZHZPWOyAi0lC9916s3K9fdHEcKfWsi0gBq6ln/R2CXvQvuPup7v7fQEUd2h4NlLr76nAmmVlA4jj3CcBDYfkxYJyZWbh+lrvvdfc1QGnYXk1tfgOY5u4HAdx9Yx1iFRFpWPbsgXXrgnKjRlUT3nzTvTsUFQXlDRtg9+5o4xERqUc1JeuXABuAF83sD2Y2jmDMerq6AevilsvCdUn3cfcDBMNtOtRQt6Y2+xL0yi80s2fMLI8/8xUROULxPdA9e0LTptHFcqQaN646XePatZGFIiJS31IOg3H3x4HHzawVcCHwXaCzmf0OeNzd/15L28kSe09zn1Trk/1zUdlmM6Dc3U8ys4uB+4HTqgVlNhmYDNCpUydKSkqSBi+Fa+fOnboupJp8uy46vPIKQ8Py1g4deCuPYk9mWLt2HB2W33rySbZuzI0PT/PtupD6oetCMqnW2WDcfRfwMPCwmbUHLgOmArUl62UEY8grdQfWp9inzMwaA22BrbXUTbW+DPjfsPw48ECK87kPuA+guLjYx44dW8tpSKEpKSlB14Ukyrvr4o03DhXbjx6dX7EnM3IkLFoEwAmtWkGOnE/eXRdSL3RdSCal81CkQ9x9q7v/3t3PSmP3BUB/M+ttZk0Jbhidk7DPHOCasHwp8IK7e7h+YjhbTG+gPzC/ljafACrjOgN4ty7nJiLSoMTPBJPPN5dWip/N5l39eheRwpHOPOuHxd0PmNkU4FmgCLjf3ZeZ2TRgobvPAWYAfzKzUoIe9Ylh3WVm9giwHDgAfNPdKwCStRke8jaC3v/vAjuBr2Tr3EREcl5Dmbax0sCB0KMHFBcHZRGRApG1ZB3A3ecCcxPW3RxXLicYVpOs7nRgejpthuu3E8wFLyIiDWXaxkpf+ELwEhEpMFlN1kVEJCLPPx8k7KWl+T1to4hIgVOyLiLSEPXr1zB61EVEClydbjAVEREREZH6o551ERHJDytXBsN73n0XTj4ZJk6MOiIRkaxTsi4i0tB89BEcfTQ0bx51JJn14oswZUpQ3r5dybqIFAQNgxERaWi++EVo1Qr69oX586OOJnMGDIiVV66MLg4RkXqkZF1EpKF55x04eBBWr4YOHaKOJnMSk3X36GIREaknStZFRBqSTz6BDRuCctOm0KtXpOFkVLdu0LJlUN6+HbZsiTYeEZF6oGRdRKQhiR8e0r8/FBVFF0ummWkojIgUHCXrIiINyTvvxMoDB0YXR7YUF8fK774bXRwiIvVEybqISEMS39vcEJN19ayLSIFRsi4i0pDE96zH90I3FOpZF5ECo2RdRKQhaejDYNSzLiIFRsm6iEhDceAAlJbGlhtiz3p8sl5aChUV0cUiIlIP9ARTEZGGYu1a2LcvKB97LBx1VKThZEXbtnDNNdCjBwweHMwn35BmvBERSaBkXUSkodi4ETp3ho8/bphDYCo9+GDUEYiI1Bsl6yIiDcUppwQPRNq2LXg4koiI5D0l6yIiDc3RRwcvERHJe7rBVERE8pd71BGIiGSVknUREckv+/fDpEkwejR06RLcZCoi0kBpGIyISEOwaRP87//CkCHBLCkNeRhMkyYwd25wQy0Es+D06RNpSCIi2aKedRGRhmD+fPjGN+C00+DCC6OOJvsGDYqVly2LLg4RkSxTsi4i0hDEJ6xDhkQXR30ZPDhWXr48ujhERLJMybqISEOwdGmsHJ/INlTx56iedRFpwLKarJvZeDNbaWalZjY1yfZmZjY73D7PzHrFbbspXL/SzM6tQ5v/bWY7s3VOIiI5qdB61uOHwahnXUQasKwl62ZWBNwNnAcMAq4ws0EJu10HbHP3fsCdwO1h3UHARGAwMB64x8yKamvTzE4C2mXrnEREclJFRdWEtdB61les0IwwItJgZbNnfTRQ6u6r3X0fMAuYkLDPBOChsPwYMM7MLFw/y933uvsaoDRsL2WbYSL/S+D7WTwnEZHcs3o1lJcH5S5doEOHaOOpDx07wjHHBOXdu2HNmmjjERHJkmwm692AdXHLZeG6pPu4+wHgE6BDDXVranMKMMfdP8pQ/CIi+aHQhsBUOuGEWHnJkujiEBHJomzOs25J1iU+ai7VPqnWJ/vnws2sK3AZMLbWoMwmA5MBOnXqRElJSW1VpMDs3LlT14VUk8vXRc+nnqJ3WC5r25bSHI0z0/q2b0+PsLz2iSdY2759vceQy9eFREfXhWRSNpP1Mjj0exSgO7A+xT5lZtYYaAtsraVusvUjgH5AaTCKhpZmVhqOha/C3e8D7gMoLi72sWPHHs65SQNWUlKCrgtJlNPXxd13Hyp2Hz+e7rkaZ6aVlcEjjwDQq7ycXhGcd05fFxIZXReSSdlM1hcA/c2sN/AhwQ2jX0zYZw5wDfAacCnwgru7mc0B/mJmdwBdgf7AfIIe92ptuvsyoEtlo2a2M1miLiLSIL35Zqw8YkR0cdS3cePg8cdh2DDo1SvqaEREsiJrybq7HzCzKcCzQBFwv7svM7NpwEJ3nwPMAP5kZqUEPeoTw7rLzOwRYDlwAPimu1cAJGszW+cgIpLz3OFLX4JFi+DttwtjJphKxx5bGE9rFZGCls2eddx9LjA3Yd3NceVygrHmyepOB6an02aSfVofTrwiInnHDG6+ufb9REQkL+kJpiIiIiIiOUrJuoiI5Lf9++Gtt+C116KOREQk45Ssi4hI/po/H1q3Dm4ynTIl6mhERDJOybqISL7atw9GjYJJk+A3vwluNi00/foF7wMED4favz/aeEREMkzJuohIvlq+HBYuhAcfhDvvDG42LTTt28NxxwXlvXuD90REpAFRsi4ikq8WL46VC2l+9UQnnRQrL1gQXRwiIlmgZF1EJF8tWhQrF3KyPmpUrKxkXUQaGCXrIiL5at68WDm+d7nQxCfr8+dHF4eISBYoWRcRyUd791YdBjN6dHSxRC3+H5W334Y9e6KLRUQkw5Ssi4jko8WLY7Og9OsHHTpEG0+U2raF4uKgXFEBb74ZbTwiIhmkZF1EJB/FD4EZMya6OHJF/CcLGrcuIg2IknURkXykZL0qjVsXkQaqcdQBiIjIYVCyXtXo0XDUUcH49fjEXUQkzylZFxHJN1u3wqpVQblpUxg2LNp4csGoUbBtGzTSB8Yi0rAoWRcRyTft28OGDUHvelkZNGsWdUTRU5IuIg2UknURkXzUuTNccEHUUYiISJapK0JEREREJEcpWRcRkYbBHR59FL71reBG0927o45IROSIaRiMiEg+WbsWPvggmP2kefOoo8ktZnDzzfDOO8Hy/PkwdmykIYmIHCn1rIuI5JM//QnOOAPatYNf/jLqaHLPaafFyv/6V3RxiIhkiJJ1EZF88tJLwde9e+HYY6ONJRedfnqsrGRdRBoAJesiIvli3z549dXY8hlnRBdLrorvWX/1VThwILpYREQyQMm6iEi+WLAA9uwJyn36QI8e0caTi3r2jL0vu3YF75mISB5Tsi4iki8qh8CAetVr8rnPxcrPPhtdHCIiGZDVZN3MxpvZSjMrNbOpSbY3M7PZ4fZ5ZtYrbttN4fqVZnZubW2a2cPh+qVmdr+ZNcnmuYmI1Lv4xFOznKR27rmx8t/+Fl0cIiIZkLVk3cyKgLuB84BBwBVmNihht+uAbe7eD7gTuD2sOwiYCAwGxgP3mFlRLW0+DAwEhgItgK9k69xEROrd9u3wyiux5fiEVKr63OegUfjnbcEC2Lo12nhERI5ANnvWRwOl7r7a3fcBs4AJCftMAB4Ky48B48zMwvWz3H2vu68BSsP2Urbp7nM9BMwHumfx3ERE6tfzz0NFRVAeORI6d442nlzWoQOMGhWUDx4M3jsRkTyVzYcidQPWxS2XAWNS7ePuB8zsE6BDuP71hLrdwnKNbYbDX64GvpMsKDObDEwG6NSpEyUlJWmfkBSGnTt36rqQaqK+LopnzKByosa1gwaxVtdojXoVF9Nr3jwA3p8zhzXHHJOV40R9XUhu0nUhmZTNZN2SrPM090m1PtknAYlt3gP8092TTrDr7vcB9wEUFxf7WI37lAQlJSXoupBEkV4X7nDFFYcWe11/Pb0+85loYskXnTsHc66fcw49e/SgZ5YOo98XkoyuC8mkbCbrZUD8vGLdgfUp9ikzs8ZAW2BrLXVTtmlm/wl0Ar6WgfhFRHLD/v3wf/8vPPMMvPVWbIiHpHb88cFLRCTPZXPM+gKgv5n1NrOmBDeMzknYZw5wTVi+FHghHHM+B5gYzhbTG+hPMA49ZZtm9hXgXOAKdz+YxfMSEalfTZvC9dfDU0/B6tVQVBR1RCIiUk+y1rMejkGfAjwLFAH3u/syM5sGLHT3OcAM4E9mVkrQoz4xrLvMzB4BlgMHgG+6ewVAsjbDQ94LvA+8Ftyjyv9z92nZOj8RkUgoURcRKSjZHAaDu88F5iasuzmuXA5clqLudGB6Om2G67N6LiIikofcg6FDjz0GV18NAwZEHZGISJ3oCaYiIrnso4+C6Qfl8Hz72zB8ONx6K8yeHXU0IiJ1pmRdRCSXnX8+9OgRJJ0bNkQdTf459dRY+bHHootDROQwKVkXEclVy5fDokWwfj38/vfQvHnUEeWff/u32Pv21lvBeyoikkeUrIuI5KqHHoqVL7gA2rWLLpZ81aZN8OlEpRkzootFROQwKFkXEclFFRXw5z/Hlq+5JvW+UrPrrouV/+d/YO/e6GIREakjJesiIrnomWeC4S8QPI3z3HOjjSefnX02HHdcUN68GZ58Mtp4RETqQMm6iEguuvvuWPlLX4ImTaKLJd8VFcGXvxxb/uMfo4tFRKSOlKyLiOSa996Dv/0tKJvBN74RbTwNwaRJwXsJ8NxzwXssIpIHlKyLiOSa//7vWPn886F37+hiaSiOOy6YGabSHXdEF4uISB0oWRcRySUffwx/+ENsecqU6GJpaG64IVZevz54uqmISI5Tsi4ikkvuvRfKy4PyyJHBzZGSGWeeCT/8YTB3/ZNPxobFiIjksMZRByAiInF+8INg9pfbbgsSSyWUmWMG06dHHYWISJ0oWRcRySXNm8PXvx7MDV5UFHU0IiISMQ2DERHJRU2aQCP9is66TZtg1aqooxARSUl/CUREonbwIOzbF3UUhaWiAn7/eyguDuaxP3gw6ohERJJSsi4iErXf/ja4mfT116OOpHCsXQvf+hZs2wavvgq//nXUEYmIJKVkXUQkSosXw403wtKlcMop8PTTUUdUGPr2halTY8s33QRLlkQXj4hICkrWRUSismULXHZZbAjMsGEwbly0MRWSH/8YTjwxKO/dCxddBJs3RxuTiEgCJesiIlHYsydIDktLg+VWrWDWrGA2GKkfTZvCww9D69bB8po1wfdk9+5o4xIRiaNkXUSkvu3eDRMmwL/+FSybwf/8T3Czo9SvgQPhz3+OLb/8Mlx4IezcGV1MIiJxlKyLiNSndevg9NPhuedi637xC7j44uhiKnQTJsAdd8SWn3su+B59+GF0MYmIhJSsi4jUl6efhlGj4I03Yut++lP43veii0kC3/1u8L2o9OabwT9Q7tHFJCKCknURkfqxeTNcfjl8/HGw3Lgx3H13cJOj5IYf/xj++MfgewNB8m4WbUwiUvCUrIuI1IeOHeHWW4Nyp07w/PNw/fXRxiTVXXcd/O1vwbSO55xTddumTcEQmYqKaGITkYKU1WTdzMab2UozKzWzqUm2NzOz2eH2eWbWK27bTeH6lWZ2bm1tmlnvsI33wjabZvPcRESq2L8fVqyAGTPgqqtg7NjqQyimTIGbb4bly+GMMyIJU9Iwbhz8/OfV1z/0UJDAH3ssfPGL8MADtFqzJvjei4hkSeNsNWxmRcDdwNlAGbDAzOa4+/K43a4Dtrl7PzObCNwOXG5mg4CJwGCgK/C8mQ0I66Rq83bgTnefZWb3hm3/rsYYDx6ETz5J74SOOqrqx6EVFfDpp+nVBTj66KrL+/fDjh3p1W3UCNq1q7pu79706zdpAm3bVl23ezfs2lV932TjM5s1q15/x470pzdr2RLatKm6bvv2YOq6dLRpE5tardLmzek/nr1duyCGeBs2wIEDSXdvtmkTlJXFVnTqFLwH8crKqr9Xqca2Hnts8D2odPBgcJNhunr0CK6BSvv3V42vJmbQq1fVdXv2wEcfpVe/SZPg+PF27IgN5ahNixbQrVvVddu2JZ/LOtn716ZN8P7F27gxaCMdRx8NxxxTdd2HH6b/s9ulS/Wf3Tlz4P33g2voo4+C1wcfwHvvVU/aFi0KnkxaqXFj+MlP0ju25BZ3uP/+oLxpE8ycCTNnMgrga1+DAQOCn5Xu3aFrV7j6aujXr2obixYFP8tNmgSvxo1j5cq/L5VfO3SAoqJY3YMHY9d94r7JvrZpU/1vVrLf+clU1o934ED6v/MbNar+O3vfPigvT69+UVEwlWm8vXuDVzqaNAl+98Tbsyf9f6qaNq0+heru3Sn/ZlSTbPrVXbvS/0SmRYuqfzMg+L2b7v0TLVvGhnJVqku+0rp11b85Bw/WbXYkXXtV11Vee02PoA/Z3bPyAk4Gno1bvgm4KWGfZ4GTw3JjYDNgiftW7peqzbDOZqBxsmOneo0MLv30Xp9+6lW88076dVu39mpefDH9+gMGVK//8MPp1z/jjOr1f/Wr9OtfcUX1+jfckH79G26oXn/ixPTr/+pX1eufcUb69f/yl+r1+/dPv/6LL1av37p1+vVXrqxa99NP06+ra696/fq89n7727i37cWgMHp0+vWnTq1+fMlPe/a4f+c77p06pfe9f/756m20aJH+tbNmTdW6mzfX7fdGeXnV+m+/nX7d9u2rx/7ss+nXHzy4ev0HHki//tlnV6//s5+lX/+aa6rXnzLlyH5uL7kk/fq/+U3s90Wlk09Ov/6jj1Y/fs+e6dd/5ZXq9Zs21bUX9bV3ww0OLHSve06dtZ51oBsQ331YBoxJtY+7HzCzT4AO4frXE+pWds8la7MDsN3dDyTZvwozmwxMBhiZbIcU/vWvf1ER1zvb4oMPqp1MKgcqKni5pKTKunaLFzM8zfq79+xhfkL9Y5YvZ1Ca9bdv387ihPrdV62iX/Ldq/n4449ZkVC/77p19Ei+ezXr1q1jVUL94zdupHOa9UtXraIsof7w7dtpl3z3apYvX87GhPqj9+yhZfLdq1m8eDHbE9adWlGR9g/PvHnz2LN+/aHlot27OS3NuqBrL8prb9Xbb7MurL9z505KSkoYtn8/R6fYv/yYY9jVpw/bhw1j+4gR7OjXDxKOL3nswgvhggtovWoVRy9cSNu336blqlW03Lix2q5vvPceO+J7xoHT9+1Le+zp66+/TvnatYeWm3zyCZ+tQ6gvvfQSHteT12rNmuBTgDTsP3CAVxKu26OXLGFYmvV37drFgoT6Xd55h4Fp1t+6bRtvJdQ/bvVq+qRZf8OGDbyTUL/fhx/SPc3673/wAWsS6g/etIlOadZ/77332Nm7NyVxbYz49FPapq5SxbJly9jUsWOVdZ8pLyfdx6UtWrSITxM+eT7d/bCvvcaffMKpadYFXXuprr11dflEPUE2k/Vkt9B7mvukWp/sWqtp/+or3e8D7gMYWVTk1T4uSeG0006r+tHMe+9VH5qSQuPWrRk7dmzCysbVP15PoWXnztXrb9kSfEyahnY9e1avv3x5cMNbMgmzH3Tu14/OifVffLH68IIUehx/PD0S68+eHcSQhn4nnEC/xPp9+wbDIdIw6MQTGZRYv1ev4P/gJPbu3UuzuGEvw8eMgZNPrrrTcccl/1gvycwRYz77WegT96O+c2dQP02nnX561Y/1Vq+uPrQlhcatWlX/3jdvDr17p1W/Ze/e1et/+mnw/qeh3eDB1euvWlV9eEClxGtv2LDq197rr8OSJWkdv8eJJ1a/9p56CuL+eapJ36FD6RvWLykpCc7lqqtgzJhgiEyXLsEwna5dYcAAmh91FM0Jeg+kATvrLPjqV4HwuhgxIvi5LCsLhrh99BEjL7wwuD7ijRgRfCS/f3/wOnAgVoYqv5M+c8opVX9PbNkC7dvH9qnl6xlnnFF1+N4xx1QfXpBCk3btqv/cVlSkXb/VMcdUr//hh2nXb9+tW/X6b7yRdv0uvXrRJbH+00+nXb9n//70TKx/3HFp1+8/aBAfJv7d79IlGC6XhsHDhgX3vMTr2DHtoSAnjhoV/I6Kd9RRaQ8d/czJJ0PPnrEVW7emfe6gay/Vtdcj1d+9NJinSFiOlJmdDNzi7ueGyzcBuPvP4/Z5NtznNTNrDGwAOgFT4/et3C+sVq1N4DZgE9Al7KGvcuxUiouLfeXKlZk4XWlADiVlInF0XUgyui4kGV0XkoyZveHuJ9W1XjZng1kA9A9naWlKcMPonIR95gDXhOVLgRc8+O9hDjAxnC2mN9AfmJ+qzbDOi2EbhG0+mcVzExERERHJuqwNgwl7uKcQ3BxaBNzv7svMbBrBAPs5wAzgT2ZWCmwlSL4J93sEWA4cAL7p7hUAydoMD/kDYJaZ3Qq8GbYtIiIiIpK3sjlmHXefC8xNWHdzXLkcuCxF3enA9HTaDNevBkYfYcgiIiIiIjlDTzAVEREREclRStZFRERERHKUknURERERkRylZF1EREREJEf9//buNVaOuozj+PcHqLywYEL7AgWtgVZFxNZUAiEGGhoDfUFfSAwEVEyjMREQvCTeEon6ArxANEG8IEGNN0SDjYiQKIqpFtsEuWpNRYNVTCUlaCQAAAZaSURBVKtWUgKi2McXO9XluOwZCGd2zu73k2zO/Gf+O/ucnCczz/nPf3Ys1iVJkqSesliXJEmSespiXZIkSeqpDB7+OZuS7AW2TzoO9c5S4C+TDkK9Y15oFPNCo5gXGuUlVbXkqb5pQR+KtAhsr6o1kw5C/ZJkm3mhucwLjWJeaBTzQqMk2fZ03uc0GEmSJKmnLNYlSZKknpr1Yv3zkw5AvWReaBTzQqOYFxrFvNAoTysvZvoGU0mSJKnPZn1kXZIkSeqtmSjWk5yWZHuSHUneO2L7c5J8s9l+e5Ll3UeprrXIi3cmuS/JXUl+mORFk4hT3ZovL4b6nZmkkviND1OuTU4keX1zvLg3yde6jlHda3EOeWGSW5Pc0ZxH1k8iTnUryTVJdiW550m2J8mnm7y5K8mr5tvn1BfrSQ4ErgROB44Bzk5yzJxuG4E9VXU0cAVwWbdRqmst8+IOYE1VHQdcD3ys2yjVtZZ5QZIlwIXA7d1GqK61yYkkK4D3ASdV1cuBizoPVJ1qeaz4IHBdVa0GzgI+022UmpBrgdPGbD8dWNG83gpcNd8Op75YB44HdlTV/VX1T+AbwIY5fTYAX2qWrwdOTZIOY1T35s2Lqrq1qh5pmluAIzqOUd1rc7wA+AiDf97+0WVwmog2OfEW4Mqq2gNQVbs6jlHda5MXBRzSLB8K/KnD+DQhVXUb8LcxXTYAX66BLcDzkhw+bp+zUKy/APjDUHtns25kn6p6HHgIOKyT6DQpbfJi2EbgpgWNSH0wb14kWQ0cWVXf6zIwTUybY8VKYGWSzUm2JBk3qqbp0CYvLgHOTbIT+D5wQTehqeeeav0xE08wHTVCPvcrcNr00XRp/TdPci6wBjh5QSNSH4zNiyQHMJgqd15XAWni2hwrDmJwSfsUBlfgfprk2Kr6+wLHpslpkxdnA9dW1SeTnAh8pcmLfQsfnnrsKdecszCyvhM4cqh9BP9/Keq/fZIcxOBy1bhLGFr82uQFSdYBHwDOqKrHOopNkzNfXiwBjgV+nOT3wAnAJm8ynWptzyHfrap/VdXvgO0MindNrzZ5sRG4DqCqfg4cDCztJDr1Wav6Y9gsFOtbgRVJXpzk2Qxu8tg0p88m4E3N8pnAj8ovoJ928+ZFM93hcwwKdeegzoaxeVFVD1XV0qpaXlXLGdzLcEZVbZtMuOpAm3PIDcBagCRLGUyLub/TKNW1NnnxAHAqQJKXMSjWd3capfpoE/DG5lthTgAeqqoHx71h6qfBVNXjSc4HbgYOBK6pqnuTfBjYVlWbgC8yuDy1g8GI+lmTi1hdaJkXHweeC3yrud/4gao6Y2JBa8G1zAvNkJY5cTPw2iT3Af8G3lNVf51c1FpoLfPiXcAXklzMYJrDeQ4ETr8kX2cwJW5pc7/Ch4BnAVTVZxncv7Ae2AE8Arx53n2aN5IkSVI/zcI0GEmSJGlRsliXJEmSespiXZIkSeopi3VJkiSppyzWJUmSpJ6yWJckSZJ6ymJdkqZYksOS/LJ5/TnJH4faP1ugz1yd5Oox25cl+cFCfLYkTZupfyiSJM2y5uE8qwCSXAI8XFWfWOCPfT/w0TEx7U7yYJKTqmrzAsciSYuaI+uSNKOSPNz8PCXJT5Jcl+Q3SS5Nck6SXyS5O8lRTb9lSb6dZGvzOmnEPpcAx1XVnU375KGR/Dua7QA3AOd09KtK0qJlsS5JAngl8A7gFcAbgJVVdTxwNXBB0+dTwBVV9Wrgdc22udYA9wy13w28vapWAa8BHm3Wb2vakqQxnAYjSQLYWlUPAiT5LXBLs/5uYG2zvA44Jsn+9xySZElV7R3az+HA7qH2ZuDyJF8FvlNVO5v1u4DnP/O/hiRNF4t1SRLAY0PL+4ba+/jfueIA4MSqepQn9yhw8P5GVV2a5EZgPbAlybqq+nXTZ9x+JEk4DUaS1N4twPn7G0lWjejzK+DooT5HVdXdVXUZg6kvL202reSJ02UkSSNYrEuS2roQWJPkriT3AW+b26EZNT906EbSi5Lck+ROBiPpNzXr1wI3dhG0JC1mqapJxyBJmiJJLgb2VtW471q/DdhQVXu6i0ySFh9H1iVJz7SreOIc+CdIsgy43EJdkubnyLokSZLUU46sS5IkST1lsS5JkiT1lMW6JEmS1FMW65IkSVJPWaxLkiRJPfUfnziEGhfou8AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Estimate Courant criterion for 5pt-FD operator\n",
    "# ----------------------------------------------\n",
    "\n",
    "dx = 0.5 # grid point distance in x-direction\n",
    "zeta = 1.0\n",
    "dt = dx / (zeta*vp0)\n",
    "op = 5 # use 5-point operator\n",
    "# RUN FD CODE WITH 5-POINT OPERATOR HERE!\n",
    "\n",
    "# PLOT SOLUTION OF 5-POINT FD OPERATOR HERE!\n",
    "\n",
    "# plot analytical solution\n",
    "Analy_seis = plt.plot(time,seis_ana,'r--',lw=3,label=\"Analytical solution\") \n",
    "\n",
    "plt.xlim(time[0], time[-1])\n",
    "plt.title('Seismogram comparison')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* Use the Taylor series to calculate finite-difference operators\n",
    "* How to derive high-order finite-difference operators\n",
    "* Derivation and application of spatial 5-point operator to approximate the spatial second derivative in the 1D acoustic wave equation"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
